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SUMMARY 

Problem 

The  problem  which  is  investigated  is  the  determination  of  the 
optimal  open  loop  controls  so  that  the  torpedo's  trajectory  minimizes 
the  cost  functional.  The  basic  cost  functional  is  the  elapsed  time 
so  that  the  problem  is  referred  to  as  the  time-optimal  torpedo  control 
problem.  For  comparison  of  different  computational  methods,  a fixed- 
time problem  of  minimizing  the  miss  distance  was  solved.  The  results 
are  restricted  to  the  torpedo's  motion  in  the  horizontal  plane,  and 
the  equations  of  motion  of  a torpedo  are  simplified  so  that  the  equa- 
tions are  linear  with  respect  to  the  control  variable.  The  computa- 
tional techniques  that  were  used  to  solve  the  problem  are  the  Epsilon 
Technique,  the  conjugate  gradient  method,  and  the  method  of  linearization. 
Results 

The  results  show  that  the  best  controls  for  the  time-optimal  tor- 
pedo control  problem  are  obtained  when  the  Epsilon  Technique  is  ini- 
tialized from  the  output  of  the  method  of  linearization  to  the  corre- 
sponding fixed-time  torpedo  control  problem.  The  basic  problem 
encountered  is  the  computer  limitation  in  the  number  of  sample  points 
and  basis  functions  able  to  be  used  on  present-day  computers  due  to 
size  limitations  and  roundoff  and  truncation  errors. 

Recommendations 

It  is  recommended  that  different  ways  of  implementing  the  com- 
putational techniques  be  explored  to  determine  how  to  increase  the 
number  of  sample  points  or  of  removing  the  need  for  more  sample  points 
by  a better  choice  of  basis  functions. 
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Chapter  One:  Introduction 

In  this  report,  the  results  of  applyinci  several  different 
coKiputational  techniques  to  the  torpedo  control  problem  are  presented. 
The  torpedo  control  problem  is  defined  to  be  the  deterr.ination  of  the 
optimal  open  loop  controls  so  that  the  torpedo's  trajectory  minimizes 
the  cost  functional.  In  the  fixed  time  torpedo  control  the  cost 
functional  is  the  distance  between  I he  torpedo  and  the  target  and  the 
final  time  is  fixed.  For  the  time-optimal  torpedo  control  problem  the 
cost  functional  is  the  final  time  and  the  taroet’s  position  is  a final 
endpoint  condition  that  the  state  equation  must  satisfy.  The  motivation 
behind  choosing  this  problem  is  an  attempt  to  answer  the  following 
question.  Given  that  the  torpedo  has  a finite  amount  of  fuel  that  is 
expended  at  a constant  rate,  will  the  effective  range  of  the  torpedo  be 
significantly  increased  if  the  time-optimal  pursuit  trajectory  is 
followed?  In  order  to  answer  the  above  question  the  time-optimal 
pursuit  trajectories  of  the  torpedo  must  be  calculated.  Hence,  we  have 
the  time-optimal  torpedo  control  problem.  Since  most  computational 
techniques  are  restricted  to  fixed  time  optimal  control  problems,  the 


fixed  time  torpedo  control  problem  was  defined  to  obtain  a comparison 
I of  different  computational  methods.  Although  this  problem  would  appear 

to  have  been  investigated  in  the  past,  the  author  has  not  been  able 
to  find  any  record  in  the  literature,  by  contacting  people  at  the 
Naval  Ocean  Systems  Center  who  have  worked  on  torpedoes  and  torpedo 
control  systems  for  over  tv/enty  years  or  by  contacting  several  professors 
at  the  Naval  Postgraduate  School. 

The  computational  techniques  used  to  solve  the  torpedo  control 
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problem  are  the  Epsilon  Technique,  the  conjugate  gradient  method  and 
the  method  of  linearization  about  a known  trajectory.  The  only  method 
which  is  directly  applicable  to  time-optimal  control  problems  is  the 
Epsilon  Technique.  A new  method  of  applying  the  Epsilon  Technique  is 
given.  In  solving  the  time-optimal  torpedo  control  problem  the  Epsilon 
Technique  is  initialized  by  using  the  results  of  the  method  of  linear- 
ization about  a known  trajectory  applied  to  the  fixed  time  torpedo 
control  problem. 

Before  proceeding  further  let  us  define  some  notational  abbrevia- 
tions that  will  be  used  throughout  this  dissertation.  The  dependence 
of  a vector  on  time  will  usually  be  implied,  that  is  x x(t),  especially 
when  being  used  as  an  argument.  Partial  derivatives  will  be  denoted  by 
a subscript,  that  is  H — 

The  Epsilon  Technique  was  first  presented  by  Balakrishnan  in 
references  [1]  - [4].  Basically,  the  Epsilon  Technique  converts  a 
dynamic  optimal  control  problem  into  a nondynamic  optimal  control 
problem  by  incorporating  the  dynamic  equation  as  a penalty  function 
added  to  the  cost  functional.  Let 

x(t)  = f(t,x,u)  (1  -1 ) 

be  the  dynamic  equation  where  the  state  vector  x(t)  is  continuous  and 


differentiable,  the  control  vector  is  contained  in  a compact  subspace  of 
and  appropiate  initial  and  final  endpoint  conditions  are  given.  We 
want  to  determine  the  controls  u(t)  that  minimize  the  cost  functional 
J(t,x,u)  = C'(t^,x(t^))  + /^f  L(t,x,u)  dt  (1.2) 

The  Epsilon  Technique  forms  the  epsilon  cost  functional 

J(t,x,u)  = J(t,x,u)  + ^ /gf  ||x-f(t,x,u)||^dt  (1.3) 
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and  minimizes  it  over  all  admissible  x and  u subject  to  the  initial  and 
final  conditions  being  met.  In  this  way  one  is  no  longer  constrained  to 
satisfy  the  dynamic  equation  (I.  I).  However,  the  minimum  of  the  epsilon 
cost  functional  must  satisfy  the  dynamic  equation  due  to  the  huge  penalty 
incurred  of  the  dynamic  equation  is  not  satisfied.  Balakrishnan  has 
shown  tliat  the  minimum  of  the  epsilon  cost  functional  converges  to  the 
minimum  of  the  original  optimal  control  problem  as  epsilon  goes  to  zero. 

The  Epsilon  Technique  has  been  applied  to  several  different  problems. 
Balakrishnan  in  reference  [4]  and  Taylor  in  reference  [5]  determined  the 
minimum  time  and  trajectory  for  a supersonic  interceptor  to  reach  a given 
altitude  and  velocity.  The  results  compared  favorably  with  the  previous 
solutions  to  this  problem  although  the  final  times  vjere  slightly  differ- 
ent. Taylor  and  Constantinides  in  reference  [6]-[7l  applied  the  Epsilon 
Technique  to  the  Earth-Mars  orbit  transfer  problem,  ''ne  final  time  that 
they  obtained  did  not  match  the  previous  results.  In  reference  [8], 

Mikami  solved  the  minimum  time  problem  of  a missile  flying  to  a fixed 
point  by  the  Epsilon  Technique.  In  reference  [9],  Mikami  solved  the 
minimum  time  problem  of  a particle  subjected  to  constant  magnitude 
thrust  which  is  controlled  by  tv;o  angles.  He  compared  the  Epsilon 
Technique  to  the  conjugate  gradient  method  using  sine  basis  functions 
and  polynomial  basis  functions  in  the  Epsilon  Technique.  Good  agreement 
between  the  three  methods  was  obtained  for  the  final  time  and  the  state 
trajectory.  However,  the  controls  were  different.  The  most  recent 
application  of  the  Epsilon  Technique  was  done  by  Hewett  in  reference 
[lO].  He  developed  a new  penalty  function  to  handle  state  constraints 
and  control  constraints  and  states  that  the  convergence  of  the  method 
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is  improved  if  the  second  order  terms  in  the  Newton-Raphson  iteration 
are  included.  This  method  was  applied  to  a missile  intercept  problem,  a 
climb  performance  problem  and  an  air  combat  maneuvering  problem. 

The  conjugate  gradient  method  was  first  presented  by  Hestenes  and 
Stiefel  in  reference  [11].  The  method  was  expanded  by  Fletcher  and  Reeves 
in  reference  [12].  Lasdon,  et  al  , expanded  the  method  to  optimal  control 
problems  in  reference  [13]  and  Pagurek  and  Woodside,  in  reference  [14], 
expanded  the  method  to  handle  optimal  control  problems  with  bounded 
control  variables.  Essentially,  the  conjugate  gradient  method  is  a 
modification  of  the  method  of  steepest  descents  that  takes  into  account 
second  order  effects  to  improve  the  convergence  properties.  Since  there 
is  no  way  of  presenting  the  method  without  going  into  specifics,  discus- 
sion of  the  method  is  deferred  to  chapter  four.  The  number  of  optimal 
control  problei'is  to  which  this  method  has  been  applied  are  too  numerous 
to  mention. 

The  method  of  linearization  about  a known  trajectory  is  a very  old 
method  of  solving  nonlinear  optimal  control  problems.  Reference  [15]  is 
one  of  many  good  books  that  discuss  the  method. 

The  arrangement  of  this  dissertation  is  as  follows.  In  chapter  two 
the  torpedo  control  problem  is  defined.  The  equations  of  motion  of  a 
torpedo  are  presented  along  with  the  cost  runctionals  for  the  fixed 
time  torpedo  control  problem  and  the  time-optimal  torpedo  control 
problem.  Then,  the  questions  of  existence  of  an  optimal  trajectory  and 
of  controllability  of  the  system  are  discussed.  In  chapter  three,  the 
Epsilon  Technique  is  presented  for  the  general  optimal  control  problem. 
The  Rayleigh-Ri tz  technique  of  expanding  the  state  vector  in  terms  of 
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known  basis  functions  and  the  Newton-Raphson  technique  for  minimizing 
the  the  epsilon  cost  functional  are  presented.  Then  the  technique  is 
applied  to  the  torpedo  control  problem.  Finally,  a discussion  is 
given  of  the  basis  functions  which  are  used  in  the  Rayl eigh-Ri tz  expan- 
sion and  a comparison  is  made  between  the  different  basis  functions. 

The  conjugate  gradient  method  is  presented  in  chapter  four.  The 
method  of  linearization  about  a known  trajectory  is  given  in  chapter 
five  along  with  a discussion  on  how  the  method  is  used  to  initialize 
the  Epsilon  Technique.  In  chapter  six  both  the  maximum  principle  and 
the  epsilon  maximum  principle  are  applied  to  the  torpedo  control 
problem.  An  explicit  form  is  obtained  for  the  optimal  controls. 

Chapter  seven  contains  the  results  of  applying  the  different  computa- 
tional methods  to  the  torpedo  control  problem.  First,  the  fixed  time 
torpedo  control  problem  is  discussed  and  then  the  time-optimal  torpedo 
control  problem.  Finally,  in  chapter  eight,  the  conclusions  are 
presented . 
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Chapter  Two:  System  Equations 

In  this  chapter  the  torpedo  control  problem  is  presented  as  a 
system  of  equations  that  are  to  be  optimized.  First,  the  equations  of 
motion  of  a torpedo  are  given  in  a state  space  formulation.  Then,  the 
cost  functionals  are  presented  for  the  two  problems  that  are  considered, 
the  fixed  time  problem  and  the  time-optimal  problem.  Finally,  a discus- 
sion of  existence  of  optimal  controls  and  of  controllability  is  given. 

The  equations  of  motion  of  a torpedo  are  derived  in  reference  [16] 
and  are  presented  in  state  space  formulation  in  reference  [17].  Since 
our  main  concern  is  with  determining  optimal  pursuit  trajectories  for 
the  torpedo  only  the  equations  that  were  used  are  presented  and  the 
reader  is  referred  to  the  references  for  their  derivation.  For  our 
purposes  the  control  variables  were  chosen  to  be  the  position  of  the 
rudders  and  elevators  of  the  torpedo.  It  is  assumed  that  they  can  change 
positions  instantaneously.  With  this  assumption  the  state  space  for  the 
torpedo  motion  equations  has  twelve  independent  variables.  The  equations 
of  motion  are  nonlinear  in  both  the  state  and  control  variables.  By 
assuming  that  the  roll  of  the  torpedo  is  always  controlled  to  zero,  the 
equations  of  motion  can  be  decoupled  into  a six  dimensional  system  for 
the  horizontal  plane  and  a six  dimensional  system  for  the  vertical 
plane.  These  systems  are  referred  to  as  the  yaw  system  and  the  pitch 
system,  respectively.  Due  to  the  lack  of  gravitational  forces,  the  yaw 
system  is  simpler  than  the  pitch  system.  Therefore  , in  order  to  solve 
the  simplest  system  first  our  investigation  is  confined  to  the  yaw 
system.  By  making  a minor  approximation  the  yaw  system  becomes  linear 
with  respect  to  the  scalar  control  variable. 
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Let  (xi,X2)  be  the  position  of  the  torpedo  in  an  inertial 
coordinate  system.  Let  (x3,xij  be  the  components  of  the  velocity  of  the 
torpedo  in  a body  fixed  coordinate  system  which  is  centered  at  the 
center  of  buoyancy  of  the  torpedo  with  the  positive  X-axis  pointing 
through  the  head  of  the  torpedo  parallel  to  its  length  and  the  positive 
Y-axis  being  90°  counterclockwise  from  the  positive  X-axis.  Let  X5  be 
the  angular  rate  of  change  of  the  torpedo  and  Xg  be  the  heading  of  the 
torpedo  in  the  inertial  coordinate  system  with  a positive  angle  measured 
1 couterclockwise  from  the  X-axis.  Then  the  equations  of  motion  of  the 

j torpedo  for  the  yaw  system  are: 

i A Ki.)  (2.1 ) 

; where. 


(x) 

= X , COS(Xf, ) - X,. 

sin(xfi) 

(2.2) 

f (a) 

' x_,  sin(xf.)  + X4 

COS(Xh) 

(2.3) 

= CjX'(  + C2X..X.,  + 

C3 

(2.4) 

f..(x) 

= Cl, X 3X1.  + Ct,X3X5 

(2.5) 

fs(x) 

= C7X3X4  + CBX3X5 

(2.6) 

ff>(x) 

= X5 

(2.7) 

9i(x) 

= 0 , i = 1,2, 3, 6 

(2.8) 

94  (x) 

= C6X5 

(2.9) 

9s(x) 

= CgX^ 

(2.10) 

The  initial  and  final  position  of  the  torpedo  are  given  and  the  control 
u is  constrained  to  be  within  the  closed  set  [-U,U].  The  values  of  the 
constants  used  in  (2.2)-(2.10)  are  given  in  Table  2.1.  If  one  wishes 
to  consider  a moving  target  then  the  state  vector  can  be  modified  in 
the  following  way.  Let  (x.j.,yj)  be  the  position  of  the  target  in  the 
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Table  2.1:  Values  of  the  Constants  Used  in  the 


Torpedo  Motion  Equations 


Constant 


Cp 

^3 

^4 

^5 

^6 

^7 

^8 


Value 

-0.005365516 

1.961174 

24.49011 

-0.06383694 

-0.08583888 

0.02323317 

-0.0009094651 

-0.2250926 

-0.01158485 
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) 


we  obtain  the 


inertial  coordinate  system, 
equations : 


Lettinq  x|=xi 


-X 


T’ 


X2=x?-yi 


x;  = fi(x)  - Xj  (2.11) 

= 1"?(x)  - Yj  (2.12) 

By  redefining  the  state  vector  to  be  ^^=(xj,  X2  ,X3 ,x,,  ,xr  ,X(, ) 
we  have  that  the  initial  position  is  now  the  difference  between  the 
initial  torpedo  position  and  the  initial  target  position  and  the  final 
position  is  the  origin.  The  two  systems  are  identical  for  stationary 
targets.  The  advantage  of  the  above  transformation  is  that  for  a moving 
target  the  final  point  is  the  origin  rather  than  a curve  in  the  inertial 
coordinate  system.  However,  since  most  of  this  thesis  only  deals  with 
stationary  targets  we  will  confine  our  discussion  to  the  unmodified  yaw 
system. 

Since  only  the  Epsilon  Technique  is  directly  applicable  to  time 
optimal  control  problems  a fixed  time  problem  was  considered  as  well  as 
the  time  optimal  control  problem.  The  cost  functional  for  the  time  optimal 
control  problem  is: 

J(t,x,u)  = t^  (2.13) 

The  cost  functional  for  the  fixed  time  problem  is 

J(t,x,u)  = 100.(xi(t^)-x,^)2+100.(x2(t^)-X2^)2  (2.14) 

For  some  cases  (2.14)  was  expanded  to  include  all  of  the  components 
of  the  state  vector.  The  problem  is  to  determine  the  control  u 
which  minimizes  the  cost  functional  J(t,x,u). 

Before  attempting  to  solve  an  optimal  control  problem  one  usually 
shows  that  an  optimal  control  exists  and  that  the  system  is  controllable. 
Existence  of  an  optimal  control  is  usually  shown  by  showing  that  the 
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system  equations  satisfy  Filipov's  condition.  For  our  case  one  would 
have  to  show  that 

[x,f(x)+g(x)u]  ^ c(l  + !lx||2)  (2.15) 

is  satisfied  for  all  admissible  x and  u.  However,  the  above  condition 
is  not  true  in  general  for  our  case  since  the  left  hand  side  contains 
cubic  terms.  Several  conditions  are  present  that  seem  to  indicate  that 
an  optimal  trajectory  exists.  If  we  assume  that  the  torpedo  varies  its 
thrust  in  order  to  keep  x^  constant  then  the  dimension  of  the  state 
vector  reduces  to  five,  the  state  equations  become  linear  and  (2.15)  is 
easily  satisfied.  Also,  for  any  constant  admissible  u there  exists  a 
stable  critical  point  x^.  Although  the  above  conditions  are  not 
sufficient  to  prove  that  an  optimal  control  exists  they  do  agree 
with  the  intuitive  feeling  that  an  optimal  control  does  exist. 

A sufficient  condition  for  controllabi 1 ity  of  a nonlinear  system 
is  that  the  following  matrix,  [B,AB, . . . ,a'^B]  have  rank  n+1  where  the 
dimension  of  the  state  vector  is  n+1,  where  B=g  and  A=f^  + g^u  are 
evaluated  at  a critical  point  of  the  system.  This  condition  holds 
if  m choose  a critical  point  such  that  u / 0.  However,  if  the 
critical  point  has  u = 0 then  the  above  condition  is  not  met. 

Since  the  above  condition  is  met  for  all  critical  points  but  one 
it  is  my  feeling  that  the  system  is  controllable. 
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Chapter  Three;  Epsilon  Technique 

The  theory  of  Balakrishnan's  Epsilon  Technique  has  been  presented 
several  times  in  the  literature,  references  [l]-[3].  Since  this  thesis 
is  mainly  concerned  with  applying  the  Epsilon  Technique  to  the  torpedo 
guidance  problem  the  presentation  here  will  be  confined  to  describing 
the  technique  with  regard  to  the  epsilon  cost  functional,  the  Rayleigh- 
Ritz  expansion  of  the  state  variables  and  the  modified  Newton-Raphson 
algorithm  used  to  determine  the  coefficients  of  the  Rayleigh-Ritz 
expansion  that  minimize  the  epsilon  cost  function.  The  presentation 
will  be  for  a general  optimal  control  problem.  Then,  the  application  of 
the  Epsilon  Technique  to  the  torpedo  control  problem  is  presented. 
Finally,  a least  squares  technique  is  presented  as  a method  for  com- 
paring different  basis  functions  for  use  in  the  Rayleigh-Ritz  expansion 
of  the  state  variables  and  the  results  of  the  comparison  for  several 
different  basis  functions  are  presented 

First,  a description  of  the  general  control  problem  is  given. 

The  problem  is  to  determine  the  control  variables  u^(t)  that  minimize 
the  cost  functional 

J(t,x,u)  = 4>(t^,x(t^))  + /pf  L(t,x,u)  dt  (3.1) 

subject  to  the  state  variables,  satisfying  the  dynamic  equation 
X = f(x,u)  (3.2) 

where, 

^(0)  - ^(t^)£6()^(t^))  = {the  set  of  admissible  final  points} 

f(x,u)  is  C'  in  x and  continuous  in  ij  and  u^  is  contained  in 
the  compact  set  U. 

The  Epsilon  Technique  reformulates  the  above  dynamic  optimal  control 
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problem  as  a nondynamic  optimization  problem  by  incorporating  the 
dynamic  equation  as  a penalty  function  added  to  the  cost  functional. 

The  epsilon  cost  functional  is  given  by, 

0^(t,x,u)  = J(t,x,u)  + /gf  11  X ' l(x,u)|12  dt/2c  (3.3) 

The  epsilon  cost  functional  is  then  minimized  over  all  admissible 
X and  u^  where  the  set  of  admissible  control  variables  has  been 
expanded  to  include  relaxed,  or  chattering  controls  and  the  set  of 
admissible  state  variables  is  the  class  of  absolutely  continuous 
functions  that  satisfy  the  initial  and  final  end  conditions.  In  the 
Epsilon  Technique  one  gives  up  satisfying  the  dynamic  equation  exactly 
in  return  for  satisfying  the  end  conditions  exactly.  The  minimization 
of  the  epsilon  cost  functional  is  performed  in  two  steps.  Let  us 
def i ne , 

I(x,u)  = I 1 X ' f(x,u)  I 1-'  (3.4) 

First,  I(x,u)  is  minimized  with  respect  to  u at  each  time  step.  Then 
the  epsilon  cost  functional  is  minimized  with  respect  to  the  state 
variables  using  the  controls  obtained  in  the  first  step.  Let  u^^  and 
x^.  be  the  controls  and  state  vector  obtained  from  the  above  two  steps. 
Balakrishnan  has  shown  in  reference  [1]  that  as  i increases  the 
sequences  converge  to  an  optimal  state  vector  and  an  optimal  control 
vector  that  is  dependent  upon  epsilon.  As  epsilon  goes  to  zero  the 
optimal  controls  and  state  vector  converge  to  the  optimal  controls 
and  state  vector  of  the  original  control  problem. 

A Rayleigh-Ritz  procedure  is  used  to  reduce  the  dimensionality 
of  the  the  system  to  finite  dimensions.  The  state  vector  is  represented 
as  a vector  in  a finite  dimensional  subspace.  Then  the  problem  is  reduced 
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to  determining  the  controls  and  the  coefficients  of  the  basis  functions 
that  minimize  the  epsilon  cost  functional.  Although  several  different 
basis  functions  were  tried  and  compared  the  following  basis  functions 
were  found  to  be  as  good  or  better  than  those  tried  and  will  be  referred 
to  as  the  standard  basis  functions.  Let 

x.(t)  = A°  + + lA.  . sin(iiit/t.)  , j = l,...,n  (3.5) 

J J J T ^ J , I T 

Then  the  end  conditions  are  satisfied  exactly  by  the  choice  of  A"  and 
A°°  since  the  basis  functions  equal  zero  at  the  end  points. 

Let  A be  the  vector'  containing  all  of  the  coefficients  which  are 
to  be  optimized.  A modified  Newton-Raphson  technique  is  used  to  determine 
the  coefficients  that  minimize  the  epsilon  cost  functional.  Let  t^. , 
i = 1,...,N  be  an  equally  spaced  partition  of  the  interval  [0,t^]  where 
we  define 


A t.^^  - t.  = t^  / ( N - 1 ) 


(3.6) 

I I 

Let  us  define 

W.  • = (x.(t.)-f  .(t.,x,u))(A/t)  i = l,..,N,  j = l,..,n  (3.7) 

I J • J ' 


and 


W^,^^J(t,x,u)'/'  ) 


Then  the  epsilon  cost  functional  is  approximated  by; 

J (t,x,u)  = W^W 

The  gradient  of  the  epsilon  cost  functional  is  given  by; 


(t,x,u)  = ( v^Wl  F 


(3.8) 

(3.9) 
-T3TTT3]- 


Neglecting  second  order  terms  the  Hessian  of  the  epsilon  cost  functional 
is  given  by; 

A ^ (t,x,u)  = (v^W)^V^W)  (3.11) 

Hence,  the  coefficients  are  updated  by  the  iteration 
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(3.12) 


In  using  the  Newton-Raphson  algorithm  described  above  it  is  not  necessary 
to  compute  the  inverse  of  the  Hessian  of  the  epsilon  cost  functional  as 
shown  in  equation  (3.12).  It  is  sufficient  to  solve  the  linear  equation 
Cy  = d (3.13) 

for  y where, 

C = A^J^(t,x,u) 

^ = ^i+1  - ^i 

Then,  the  new  coefficients  are  given  by; 

l\..  = A.  + y (3.14) 

Since  C is  very  ill-conditioned  we  follow  the  suggestion  given  by 
Luenberger  in  reference  [18]  and  add  a small  constant  to  the  diagonal 
elements  of  C.  Even  with  the  addition  of  0.01  to  the  diagonal  elements  of 
the  Hessian  matrix  C the  condition  number  of  the  Hessian  matrix  was  seven 
or  eight.  The  vector  y is  computed  using  a double  precision  routine 
that  decomposes  C into  a lower  triangular  matrix  and  an  upper 
triangular  matrix  Cy  such  that  C = CyCy.  Then  y is  found  by  two  successive 
back  substitutions. 

For  the  torpedo  control  problem  the  dynamic  equation  has  the  form 
^(t)  = f(x)  + g^(x)u  (3.15) 

The  cost  functional  for  the  fixed  time  torpedo  control  problem  is 
given  by; 

vl(t,x,u)  = 100(x^ (t^)-x^y)^  + 100(x2(ty)-X2^)^  (3.16) 

and  for  the  time-optimal  torpedo  guidance  problem  the  cost  functional 
is  given  by; 
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(3.17) 


J(t,x,u)  = 

Hence  the  epsilon  cost  functional  is  given  by 

J (t,x,u)=J(t,x,u)  + - /^f  I lx-f(x^)-£(x)u|  p dt  (3.18) 

Since  the  epsilon  cost  functional  is  linear  in  the  control  variable  u 
an  analytic  expression  exists  for  the  control  variable  that  minimizes 
I(x,u)  = 1 |x  -f(x)-g(x)u| |2  (3.19) 

Taking  the  partial  derivative  of  I(^,u)  with  respect  to  u and  equating  to 
zero  the  following  expression  is  obtained  for  the  optimal  control  variable 

(x4-f4(x))g4(0  + (x5  -f5(x))g5(x) 

Uo  = (3.20) 

94(x)2  + 95(x)2 

Since  the  second  derivative  of  I(x,u)  is  equal  to  the  denominator  of 
the  above  expression  which  is  always  positive  Uo  is  the  minimum  of 
I(x,u). 

In  any  application  of  the  Epsilon  Technique  it  is  important  to  use 
basis  functions  that  can  adequately  represent  the  state  vector  over  the 
entire  time  interval. In  order  to  determine  whether  the  basis  functions 
given  in  (3.5)  are  adequate  to  represent  the  trajectory  of  a torpedo 
and  to  obtain  a comparison  of  several  different  basis  functions  the  method 
of  least  squares  approximation  was  used  to  obtain  the  best  coefficients 
for  the  basis  functions  to  approximate  a known  trajectory.  Then  the 
trajectory  obtained  by  using  these  coefficients  can  be  compared  with  the 
actual  trajectory  to  see  if  the  approximation  of  the  trajectory  is 
adequate.  In  this  way  a comparison  can  be  made  between  diffe>'ent  basis 
functions.  Now,  the  method  of  least  squares  approximation  is  presented. 
Using  the  same  partition  of  the  time  interval  given  before  let  D be  a 
matrix  with  the  i,j  ^element  being  the  value  of  the  j ^ state  variable 
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at  the  i ^ time,  B be  a matrix  with  the  i,k  element  being  the  value  of 
the  k ^ basis  function  at  the  i ^ time  step  and  A be  the  matrix  contain- 
ing the  coefficients  of  the  basis  functions  with  the  i , j tji  element  being 
the  j ^ coefficient  for  the  i state  variable.  The  coefficients  A are 

obtained  from  the  equation 

A^  = (3.21) 

In  Table  3.1  a list  of  the  basis  functions  that  were  tried  is 
given.  These  functions  were  first  tried  with  the  epsilon  technique 
for  the  general  minimum  time  problem.  The  basis  functions  numbered 
3 through  6 were  rejected  because  the  basis  functions  approach  zero  as 
time  increases  and  the  trajectory  of  the  torpedo  converges  to  the  straight 
line  connecting  the  initial  and  final  points.  The  polynomial  functions 
shown  as  number  seven  didn't  perform  well  at  all.  Therefore,  only  the 
standard  basis  functions  and  the  basis  functions  which  were  proposed  by 
Taylor  in  reference  [19]  and  are  numbered  two  in  the  table  are  compared. 
Taylor  modifies  the  standard  basis  functions  by  approximating  the  initial 
and  final  conditions  with  a cubic  instead  of  a straight  line  and  by 
multiplying  each  basis  function  by  the  term  sin(7it/t^).  The  latter 
modification  causes  both  the  state  variables  and  their  first  derivatives 
to  be  equal  to  zero  at  the  initial  and  final  times. In  order  to  understand 
the  effect  of  each  of  the  above  modifications  the  case  of  the  standard 
basis  functions  with  a cubic  approximation  to  the  initial  and  final  end 
point  conditions  was  considered.  The  three  cases  that  are  considered  for 
comparison  are;  a maximun  turn  for  0.6  seconds,  a maximum  turn  for  0.3 
seconds  followed  by  zero  rudder  deflection  for  0.3  seconds  and  a maximum 
turn  for  0.3  seconds  followed  by  9.7  seconds  of  zero  rudder  deflection. 


I 
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The  values  of  the  optimal  coefficients  for  the  different  basis  functions 
are  presented  in  Tables  3.2  - 3.8.  The  values  of  the  comparisons  of  the 
basis  functions  with  the  actual  trajectories  are  presented  in  Tables  3.9  - 
3.11.  The  results  show  that  the  position  of  the  torpedo  is  represented 
adequately  by  all  the  basis  functions,  however  the  velocity  components  of 
the  state  vector  are  approximated  better  by  the  standard  basis  functions. 
There  is  no  significant  difference  between  the  case  where  the  initial 
trajectory  is  a straight  line  or  a cubic.  Tlierefore,  the  standard  basis 
functions  with  a straight  line  initial  approximation  between  the  endpoints 
was  used  in  applying  the  tpsilon  Technique. 

In  applying  the  tpsilon  Technique  the  value  of  epsilon  is  usually 
decreased  until  no  improvement  is  sliown.  For  the  torpedo  control  problem 
the  results  were  virtually  identical  for  epsilon  having  a value  of  either 
d.Ol  or  U.OOl.  The  results  started  to  deteriorate  when  a value  of  0.0001 
was  used  for  epsilon  as  a result  of  computational  problems.  Therefore, 
the  value  of  epsilon  that  v;as  used  in  0.001  although  occassional ly 
epsilon  has  the  value  O.Ul. 

There  are  several  computational  problems  in  applying  the  Epsilon 
Technique  to  the  torpedo  control  problem.  The  epsilon  cost  functional 
is  approximated  in  (3.9)  by  evaluating  it  at  a finite  number  of  sample 
points  and  treating  the  epsilon  cost  functional  as  a constant  between  the 
sample  points.  At  each  of  these  sample  points  the  gradient  of  the 
epsilon  cost  functional  with  respect  to  A is  evaluated.  Therefore,  a 
large  amount  of  computer  storage  is  needed  in  order  to  apply  the  Epsilon 
Technique  to  the  torpedo  control  problem.  For  twelve  basis  functions 
and  fifty  sample  points  over  100,000  words  of  computer  storage  are  needed. 
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Special  techniques  had  to  be  used  to  obtain  that  much  storage  on  an  Univac 
1110  computer  system.  In  applying  tne  Epsilon  Technique  to  the  torpedo 
control  problem  it  was  determined  that  50  sample  points  were  needed  for  a 
trajectory  of  0.6  seconds  duration  and  even  then  the  solution  was  not  exact. 
On  tiiis  basis  500  sample  points  would  be  needed  for  a trajectory  of  6.0 
seconds  duration  which  would  require  more  computer  storage  than  is  presently 
available  in  one  computer  system.  Present  day  computers  are  unable  to 
handle  that  amount  of  computer  storage  without  elaborate,  time-consuming 
techniques.  A second  problem  is  that  tfie  number  of  basis  functions  is 
limited.  For  12  basis  functions  with  a free  final  endpoint  the  dimension 
the  system  of  equations  that  has  to  be  solved,  (3.13),  is  equal  to  77. 

The  number  of  steps  needed  to  solve  the  matrix  equation  is  proportional 
to  the  cube  of  its  dimension.  The  more  steps  that  are  used  increases  the 
truncation  and  roundoff  errors  in  the  computational  algorithm.  Several 
cases  were  tried  using  14  basis  functions  but  there  was  no  improvement  in 
the  results.  Another  problem  is  the  slow  convergence  for  the  cases  where 
the  initial  point  is  fixed.  This  problem  is  typical  of  penalty  function 
techniques  and  the,  see  [18],  Epsilon  Technique  is  of  the  penalty  function 
genre.  In  order  to  try  and  increase  the  speed  of  convergence  the  method 
of  steepest  descents  was  used  for  the  first  nine  iterations  and  then  the 
Newton-Kaphson  algorithm  was  used.  Another  method  was  to  apply  Aiken's 
method  of  increasing  the  speed  of  convergence  of  a convergent  sequence. 

Both  of  these  methods  didn't  give  any  noticeable  improvement  in  the  results. 
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No. 


Basis  Functions 


r 


1 sin(i7it/t^) 

2 sin(iiTt/t^)  sin(7rt/t^) 

3 sin(i~t/t^) 

4 sin(iiTt/t^)  if  i is  odd 

e cos(ii't/t^)  if  i is  even 

5 sin(Ut) 

6 e sin(Ut)  + e cos(Ut) 

7 (t/t^)\t^-t) 

Table  3.1:  l.ist  of  the  basis  functions  that  were  considered 
for  use  with  the  Epsilon  Technique. 
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Basis  functions  are  type  No.  1 with  straight  initial  approxination • Trajector'y 
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TABLE  3.11:  Comparison  of  values  for  0.3  sec.  curve,  9.7  sec.  straight  trajectory' 
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Chapter  Four:  Conjugate  Gradient  Method 
The  conjugate  gradient  method  was  developed  by  Hestenes  and  Stiefel, 
reference  [11],  and  expanded  by  Fletcher  and  Reeves,  reference  [12],  to 
be  a rapidly  convergent  finite  dimensional  unconstrained  minimization 
technique.  In  reference  [13],  Lasdon,  fitter  and  Warren  expanded  the 
conjugate  gradient  method  to  optimal  control  problems.  Pagurek  and 
Woodside,  reference  [14],  further  extended  the  method  to  handle  optimal 
control  problems  with  bounded  control  variables.  In  their  paper  two 
separate  algoritltns  were  presented.  The  first  is  a direct  extension  of 
Lasdon's  algorithm.  The  second  algorithm  uses  second  order  terms  and 
gave  improved  convergence  in  several  sample  cases.  Therefore,  the  second 
method  was  used  for  the  torpedo  control  problem.  The  discussion  here 
will  be  confined  to  presenting  the  steps  of  the  algorithm  for  the 
general  optimal  control  system  that  has  the  form  of  the  torpedo  control 
problem.  Let  the  dynamic  equation  be  given  by 

x(t)  = f(x,u)  (4.1 ) 


where, 

21(0)  = Xo  , ut  [-U,U] 

The  cost  functional  has  the  form 
J(t,x,u)  = i})(t^,x(t^)  i 

The  Hamiltonian  of  the  above  system  is  given  by 
H(t,x,£,u)  = £^f(x,u) 

where  £ is  the  costate  variable  defined  by 
p(t)  = - 8H(t,x,£,u) 


Txltr 


, £(t^)=aJ(t,x,u) 


9X 


■f 


(4.2) 


(4.3) 


(4.4) 


Let  z(t)  and  c(t)  be  n-dimensional  vectors  and  s(t),w(t)  and  e(t)  be 
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scalars  and  let  us  denote  partial  derivatives  by  subscripts,  that  is 
means  the  partial  derivative  of  £ with  respect  to  Then,  the 
conjugate  gradient  method  is  given  in  the  following  six  steps. 


Step  1 : 

Determine  an  initial  control  u(t). 

Step  2: 

Solve  the  following  equations  forward  in  time; 

X = f(x,u)  , x(0)  = Xo 

(4.1) 

z=fz+fws  ,z(0)=0 

(4.5) 

Step  3; 

Solve  the  following  equations  backward  in  time; 

i , £(t^)  = J^(t^,x,u) 

(4.4) 

-flu  - , c(tp=J^^(t^,x,u)z(tf) 

(4.6) 

and  at  the  same  time  evaluate; 

(4.7) 

d.  = H z + H s.  + f^  c 
1 UX-  UU  1 -u  — 

(4.8) 

I^s  = .-of  we.d.  dt 

(4.9) 

Iss  = ws.d.  dt 

(4.10) 

Step  4: 

Determine  a new  value  of  s by 

^i+1  = "i  " 6i^ 

(4.11) 

where , 

6i  = ■^vs'^^ss  ' So  = 6o  ,eo  = 0 

(4.12) 

Step  5; 

Determine  a new  control  u from  the  equation 

u = u - as^. 

(4.13) 

where  a is  determined  by  a one  dimensional  search  to  be  the 
value  that  minimizes  J(t,x,u).  During  the  search  u is  truncated  to  be 
within  the  1 imits  of  U. 

Step  6:  After  defining  a new  u,  w is  defined  by 


w 


Gifu 
1 if  u 


= U 
< U 
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Then  set  i equal  to  i+1  and  go  to  Step  2. 

The  method  is  terminated  when  Step  5 yields  no  improvement  in  the  cost 
functional  or  when  a preset  number  of  iterations  have  occurred. 

For  the  torpedo  control  problem  the  dynamic  equation  has  the  form 
x(t)  = f(x)  + 3^(x)u  (4.14) 

Since  the  conjugate  gradient  method  is  not  designed  for  time-optimal 
control  problems  it  is  only  applied  to  the  fixed  time  torpedo  control 
problem.  Therefore,  the  cost  functional  is  given  by 

J(t,x,u)  = 100(x^(t^)-x^^)2  + 100(x2(t^)  - X2^)2  (4.15) 

Since  the  dynamic  equation  is  linear  in  u,  is  equal  to  zero. 

The  integrations  in  Step  2 and  in  Step  3 were  originally  performed 
using  a four  point  predictor-corrector  method  which  is  known  as  the 
Improved  Adams  or  Moulton's  method.  Several  of  the  cases  were  rerun 
using  the  simple  one  point  Euler  method.  Due  to  the  small  step  size  there 
was  no  noticeable  difference  in  the  results.  Therefore,  in  order  to 
save  computer  time  the  simple  Euler  method  was  used  for  integrating 
the  dynamic  equations.  The  one  dimensional  search  in  Step  5 was  done  by 
evaluating  the  cost  functional  at  20  equally  spaced  points  and  then 
fitting  a quadratic  about  the  minimum  of  the  20  points.  If  the  minimum 
occurred  on  the  first  point  and  the  value  of  the  cost  functional  was 
not  reduced,  a new  search  was  performed  in  the  interval  between  the 
old  point  and  the  first  point  of  the  old  search.  If  after  reducing 
the  region  of  search  several  times  no  improvement  is  made  the  algorithm 
is  assumed  to  have  converged. 
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Chaptvir  Five:  flethod  of  Linearization 
The  method  of  linearization  about  a known  trajectory  has  always  been 
one  of  the  most  frequently  used  methods  to  solve  nonlinear  optimal 
control  problems.  One  advantage  of  this  method  is  that  a feedback 
control  is  obtained.  This  method  is  only  applied  to  the  fixed  time 
torpedo  guidance  problem. 

The  trajectory  >^'(t)  about  which  the  trajectory  is  linearized  is 
a straight  line  between  the  initial  and  final  end  points.  Let  us  define 
^(t)  = )^(t)  - x'(t).  Then,  for  the  torpedo  motion  equations  which  were 
presented  previously,  the  dynamic  equation  for  the  linearized  state 
vector  is; 

z^(t)  = Cz(t)  + Du(t)  (5.1) 

where , 


z(0)  = lo  , z_(t^)  = 0 , C = 3f(_xJ  + 3g(x)u 

9x  x'(t)  !x'(t) 

D = 0(l)l,.(t) 

Although  in  general  C and  0 are  time  dependent,  for  the  torpedo 
control  problem  that  we  are  considering  with  the  initial  trajectory 
being  a straight  line  C and  D are  independent  of  time.  C and  D are 
given  by; 


0 

0 

COS(Xg) 

-sin(x^) 

0 

-X2Sin(xg) 

0 

0 

sin(xg) 

COS(Xg) 

0 

x^cos(x^) 

0 

0 

2C1X3 

0 

0 

0 

(5.2) 

0 

0 

0 

^4^3 

^5^3 

0 

0 

0 

0 

*^8^3 

0 

_0 

0 

0 

0 

1 

0 _ 
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(5.3) 


= ( 0 0 0 ® ^ 

In  order  to  insure  that  the  constraints  on  the  control  u are  satisfied 
a penalty  function  of  the  integral  of  u is  added  to  the  cost  functional 
and  the  amount  of  the  penalty  is  increased  until  the  constraints  are 

satisfied.  The  cost  functional  is  given  by; 

T ^f 

J(t,z,u)  = z Ez  +A/g  u(t)^  dt  (5.4) 

where  E has  non  zero  values  only  on  the  diagonal.  E-^  is  egual  to 
100  if  x.j(t^)  is  fixed.  Otherwise,  E^^  is  equal  to  zero.  The  Hamiltonian 
of  the  above  system  is; 

H(t,z,p,u)  = Po>u^  + p^Cz  + n^Du  (5.5) 

Letting  po  = -1/2  and  setting  the  derivative  of  the  Hamiltonian  with 
respect  to  u equal  to  zero  vie  obtain  the  following  relationship  for  u; 

u*(t)  = dV  (5.6) 

Hence,  the  optimal  state  and  costate  vectors  satisfy  the  equations; 


^(t)  = Cz(t)  + I{)D^p(t)  ( ") 

p*(t)  = -Cp(t)  (5.8) 

Let  us  assume  that  p(t)  = R(t)z(t) . Then  we  have 

p(t)  = R(t)z(t)  + R(t)z(t) 

Using  (5.7)  and  (5.8)  in  (5.9)  we  obtain 

-C^R(t)z(t)  = R(t)z(t)  + R(t)Cz(t)  + |R(t)DD^R(t)z(t)  (5.10) 

Since  z(t)  is  arbitrary  we  must  have; 

R(t)  = -C^R(t)  - R(t)C  - jR(t)DD^R(t)  (5.11) 

u*(t)  = jT)^R(t)z*(t)  (5.12) 

z*(t)  = (C  + JoD^R(t))z(t)  (5.13) 


Therefore,  the  optimal  trajectories  and  controls  are  obtained  by 
solving  (5.11)  - (5.13)  subject  to  the  endpoint  conditions  on  ^(t) 
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and  R(t^)  = -E. 

The  method  of  linearization  about  a known  trajectory  is  used  as  a 
way  of  initializing  the  trajectory  for  the  Epsilon  Technique.  First, 
the  fixed  time  torpedo  control  problem  is  solved  using  the  method  of 
1 inearization  for  an  approximate  final  time  that  is  less  than  the  actual 
final  time.  Then,  the  method  of  least  squares  approximation  is  used  to 
obtain  the  optimal  coefficients  of  the  Rayl eigh-Ri tz  expansion  of  the 
state  vector  using  the  standard  basis  functions.  The  method  of  least 
squares  approximation  is  given  in  chapter  3.  Finally,  the  minimum  time 
torpedo  control  problem  is  solved  by  the  Epsilon  Technique  using  the 
loetficients  obtained  above  as  a start  off  point. 


•s 


Chapter  Six:  Application  of  the  Maximum  Principle 
In  this  chapter  the  formalism  of  ttie  maximum  principle  is  applied 
to  the  torpedo  control  problem  to  obtain  the  expectea  form  of  the 
optimal  controls.  Then  the  Epsilon  Technique  formulation  of  the 
torpedo  control  problem  is  given  and  the  Epsilon  maximum  principle  is 
applied.  Since  the  only  differences  between  the  time  optimal  torpedo 
control  problem  and  the  fixed  time  torpedo  control  problem  is  the 
addition  of  the  constant  Po  in  the  Hamiltonian  of  the  time  optimal 
problem  and  the  way  of  initializing  the  costate  vector  p(t),  we  will 
confine  our  discussion  to  the  time  optimal  torpedo  control  problem. 
The  dynamic  equation  for  the  torpedo  control  problem  is 


xU)  = f (x)  + cj(x)u  (6.1  ) 

The  cost  functional  is 

J(t,x,u)=t^  (6.2) 

The  Hamiltonian  of  the  system  is  aiven  by 

H(t,x,p,u)  = p(t)^f(x)  + p(t)^g(x)u  + p,  (6.3) 

The  costate  vector  p(t)  satisfies  the  equation 
p(t)  = - ''<(t,X,p,u) 

■X 


Now,  it  is  a necessary  condition  of  the  maximum  principle  that  the 
Hamiltonian  is  maximized  for  the  optimal  state  vector  x*(t)  and  the 
optimal  control  u*(t).  From  (6.3)  we  see  that  the  optimal  control  must 
have  the  form 

u*(t)  = sgn(p*(t)^g(t))  (6.5) 

provided  that  p*(t)g(t)  is  not  equal  to  zero  for  any  interval  of 
time.  If  H^(t,x*,p*,u*)  = p*(t)^g(t)  = 0 over  a non  zero  interval  of 
time  then  we  have  singular  controls.  For  this  to  be  the  case  we  miist 
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have  H^^(t,x*,p*,u*)  = 0.  Since  (6.1)  is  autonomous  we  must  have 
H(t,x*,p*,u*)  = 0 for  the  time  optimal  case  by  the  maximum  principle. 
Using  the  above  equations  we  can  obtain  expressions  for  P^.P^  and  p^  in 
terms  of  p.|,p2,p^  and  x*.  If  we  take  the  above  expression  for  p^  and 
differentiate  it  with  respect  to  time  and  try  to  equate  it  with  the 
expression  for  p^  we  obtain  four  relationships  that  must  hold. These 
relationships  are  consistent  only  if  u = x^  = x^  = 0.  Therefore,  the 
only  possible  singular  control  is  u = 0 with  the  torpedo  moving  in  a 
straight  line.  Therefore,  the  expected  optimal  controls  are  *U  and  0. 

Wow,  let  us  consider  the  epsilon  maximum  principle.  Define, 


p^  (t)  = ^x'  (t)  - f(x'  ) - g(x'  )u'  ) (6.6) 

where  \ and  u are  the  state  vector  and  controls  that  minimize  the 
epsilon  cost  functional 

J (t,x,u)  = J(t,x,u)  + ^ .'pf  l(x  - f(x)  - 9(x)u(l^  dt  (6.7) 

Then,  by  the  epsilon  maximum  principle, 

1 im  p (t)  = p(t)  (6.8) 

, -0 

Define, 

H (t,x'  ,p‘  ,u  ) = p (t)^^(  f(x'  ) + q(x'  )u’  ) (6.9) 


The  epsilon  maximum  principle  states  that  the  epsilon  Hamiltonian 
should  be  a maximum  at  the  optimal  control  and  state  vector.  From 
(6.9)  it  is  clear  that  we  should  have 

u*  = Usgn(p‘ (t)^g(x' ) ) (6.10) 

Taking  the  derivative  of  the  epsilon  Hamiltonian  with  respect  to 
the  control  we  have; 

H[j(t,x’  ,p' ,u)  = P^g^(x' ) + P595(x  ) = 0 (6.11) 
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(6.12) 


Usinq  (6.6)  in  (6.11)  and  qatherinq  tom'",  wo  obtain 

u'  = (^4  - " (^5  - *5^% 

94  ^ 95 

The  above  expression  for  u is  exactly  wrtai  was  obtained  in  chapter 
3 by  minimizing  1 (x,u) . 

Since  the  maximum  principle  gives  the  fotm  of  the  0[)timal  controls 
to  be  either  '11  or  0,  the  Epsilon  Technique  and  the  conjugate  (jradient 
method  were  t»'ied  with  the  controls  restricted  to  be  those  given  by  the 
maximum  principle.  The  results  of  limiting  the  control  in  this  fashion 
will  be  given  in  the  next  chapter. 
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Chapter  Seven:  Results 

The  results  of  applying  the  Epsilon  Technique,  the  conjugate  gradient 
method  and  the  method  of  linearization  about  a known  trajectory  to  the 
torpedo  control  problem  are  presented  in  this  chapter.  All  of  the  above 
methods  were  applied  to  the  fixed  time  torpedo  coritrol  problem.  Only  the 
Epsilon  Technique  is  applied  to  the  time-optimal  torpedo  control 
problem.  There  are  two  ways  of  using  the  Epsilon  Technique;  one  is  the 
standard  way  of  assuming  a straight  line  between  the  endpoints  as  the 
initial  trajectory  of  the  state  variables  and  the  second  way  is  to  first 
solve  a fixed  time  torpedo  control  problem  to  obtain  an  initial  estimate 
for  the  trajectory  of  the  state  variables.  The  arrangement  of  this  chapter 
is  as  follows.  First,  definitions  will  bo  given  of  several  terms  that 
will  be  used  throughout  the  text.  Then  the  fixed  time  torpedo  control 
problem  will  be  discussed.  The  last  part  of  tiie  chapter  is  devoted  to 
the  time-optimal  torpedo  control  problem. 

In  the  torpedo  control  problem  the  state  vector  is  composed  of  six 
variables;  the  position  of  the  torpedo  in  the  inertial  coordinate  system, 
the  velocity  of  the  torpedo  in  the  body  fixed  coordinate  system  and  the 
angular  velocity  and  heading  of  the  torpedo.  Only  the  first  two  components 
of  the  state  vector,  those  referring  to  the  position  of  the  torpedo,  are 
always  held  fixed  at  the  initial  and  final  endpoints.  By  the  term  fixed 
endpoint  we  mean  that  the  last  four  components  of  the  state  vector  are 
held  constant  and  by  the  term  free  endpoint  we  mean  that  the  last  four 
components  of  the  state  vector  are  allowed  to  vary  in  computing  the 
optimal  trajectories.  The  term  x-trajectory  refers  to  the  trajectory  of 
the  state  vector  generated  by  the  optimization  method.  The  term 
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u-trajectory  refers  to  the  trajectory  of  the  state  variables  generated 
from  the  controls  (jiven  by  the  optimization  method.  A curve  is  defined  to 
be  the  trajectory  generated  by  having  the  rudders  of  the  torpedo  at  tlieir 
maximum  value  and  a straight  trajectory  is  the  trajectory  generated  by 
having  the  rudders  with  no  deflection  at  all.  The  optimization  methods 
that  are  used  are  referred  to  by  the  followitr.i  abbreviations: 

SE  - the  standard  Epsilon  Technique 

1 

1 

ME  - the  Epsilon  Technique  with  the  control  restricted  to  be  either  i 

♦U  or  0. 

CG  - the  conjugate  gradient  method 

TiCG  - the  CG  method  with  the  contols  the  same  as  the  ME  method. 

i 

Lin  - the  method  of  linearization  about  a known  trajectory  between 
the  endpoints. 

LSE  - the  method  of  using  the  output  of  the  LIN  method  to  initialize 
the  trajectory  of  tne  state  vector  for  the  SE  method. 

From  the  results  of  appl/ino  the  maximum  principle  to  the  torpedo  control 
problem  one  has  that  the  optinial  controls  should  be  either  ‘U  or  0. 

Therefore,  several  cases  wore  chosen  for  comparison  where  the  final 
endpoint  was  obtained  by  a curve  or  by  a curve  follov^ed  by  a straight 
trajectory.  In  this  way  one  can  obtain  an  estimate  of  how  well  the 
various  optimization  algorithms  work  in  solving  the  torpedo  control 
problem.  The  cases  that  were  chosen  are: 

case  A - an  0.3  second  curve. 

case  B - an  0.6  second  curve. 

case  C - an  0.3  second  curve  followed  by  an  0.3  second  straight 

trajectory. 
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case  D - an  3.0  second  curve. 

case  E - an  0.5  second  curve  followed  by  an  2.5  second  straight 

trajectory. 

Now  the  results  of  the  fixed  time  torpedo  control  problem  are 
presented.  Each  case  will  be  discussed  separately.  Then  the  results  will 
be  summarized  for  general  conclusions.  Table  7.1  and  Figures  7. 1-7. 8 
contain  the  results  for  case  A.  As  could  be  expected,  the  ME  method  gives 

an  x-trajectory  and  a u-trajectory  which  are  exactly  the  same  as  the 

original  trajectory  from  which  the  endpoints  were  obtained.  The  SE  method 
gives  a very  good  approximation  but  the  controls  take  too  long  to  reach 
the  maximum  and  hence  the  u-trajectory  does  not  reach  the  desired 
final  endpoint.  The  CG  method  gives  a good  trajectory  however  the 
control  changes  to  +(J  from  -U  at  the  end  of  the  run  and  this  seems 
to  only  affect  the  heading  at  the  end  of  the  run  and  not  the  position. 
Several  cases  were  run  with  the  LIN  method.  For  all  of  the  cases  the 
controls  decreased  to  zero  and  the  resulting  trajectories  didn't 
develop  the  proper  curvature.  Hence,  the  LIN  method  gives  the  worst 
results  for  this  case. 

In  addition  to  comparing  the  different  optimization  algorithms, 
case  B is  used  to  compare  different  ways  of  implementing  the  algorithms. 
The  results  are  presented  in  Table  7.2  and  Figures  7.9  - 7.22.  There 
are  three  separate  cases  that  are  considered;  the  torpedo  initially 
heading  straight,  the  torpedo  initially  in  a curve  and  the  torpedo 
initially  heading  straight  with  the  aimpoint  being  the  head  of  the 
torpedo  instead  of  the  center  of  buoyancy  of  the  torpedo.  First,  we 
discuss  the  case  of  the  torpedo  initially  heading  straight. In  the  SE 
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method  the  buildup  of  the  angular  velocity  is  gradual  and  since  the  '' 

controls  depend  upon  the  angular  velocity  the  optimal  controls  are  not  i] 

reached  immediately.  Therefore,  even  though  the  x-trajectory  is  almost  | 

. i 

exact  the  u-trajectory  is  slightly  flatter  and  doesn't  reach  the  final  j- 

endpoint  or  the  desired  heading.  If  the  final  endpoint  is  fixed  the 

controls  reach  the  maximum  value  sooner  and  both  the  x-trajectory  and 

the  u-trajectory  are  closer  to  the  optimal  trajectory.  The  CG  method 

has  good  results  for  the  position  of  the  torpedo  but  the  controls  change 

sign  abruptly  at  the  end  of  the  run  and  the  heading  is  slightly  off. 

If  the  final  endpoint  is  fixed  the  control  stays  the  same  for  the  entire 

run  and  the  result  is  the  exact  trajectory.  The  CG  method  contains 

several  differential  equations  that  are  solved  numerically.  Case  B was 

used  to  compare  solving  these  equations  by  a four  point  predictor- 

corrector  routine  known  as  Moulton's  method  with  solving  these  equations 

by  the  simple  Euler's  method.  As  can  be  seen  in  Table  7.2  and  Figures 

7.15  and  7.16  the  results  are  nearly  identical.  Therefore,  the  simple 

Euler's  method  was  used  in  all  further  runs  with  the  CG  method  to  save 

computer  time.  The  ME  method  gives  the  exact  trajectory.  The  LIN  method 

gives  the  poorest  results  due  to  the  control  decreasing  in  absolute  ! 

value  to  zero.  For  the  cases  where  the  torpedo  was  initially  in  a curve 

all  of  the  methods  tried  gave  very  good  results.  This  case  was  used  to 

see  if  a noticeable  improvement  could  be  seen  in  the  way  we  applied  the 

CG  method  over  the  simpler  way  of  applying  the  CG  method.  No  improvement 

was  noticed  but  since  Pagurek  and  Woodside  showed  some  improvement  by 

using  the  more  complicated  method  we  continued  to  use  the  second  order  CG 

method.  When  we  changed  the  aimpoint  of  the  torpedo  to  be  its  head 
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i'-'tead  of  its  center  of  buoyancy  we  found  that  the  SE  method  gave  almost 
exactly  the  same  results  as  the  previous  case  where  the  aimpoint  is  the 
torpedo's  center  of  buoyancy.  However,  the  CG  method  gives  the  exact 
trajectory  and  controls.  This  seems  to  indicate  that  the  trajectory  of 
the  state  variables  are  unique  only  when  the  aimpoint  of  the  torpedo  is 
its  head  and  not  its  center  of  buoyancy. 

For  case  C the  results  are  given  in  Table  7.3  and  Figures  7.23  - 7.32. 
The  CG  method  gives  a good  estimate  of  the  final  position  but  the  final 
heading  is  slightly  larger  than  desired.  The  controls  are  always  greater 
in  absolute  value  than  0.15,  that  is  the  trajectory  is  more  of  a 
continuous  curve  than  a curve  followed  by  a straight  trajectory.  The 
MCG  method  gives  the  same  results  as  the  CG  method  gave  in  case  B.  The 
SE  method  gives  a good  approximation  of  the  x-trajectory  but  the  controls 
are  too  small  and  the  u-trajectory  is  not  close  to  the  x-trajectory. 

The  ME  method  gives  an  x-trajectory  that  is  close  to  the  desired  final 
position  of  the  torpedo  but  the  trajectory  is  even  more  curved  than  the 
trajectory  from  the  CG  method.  The  u-trajectory  is  the  same  as  the  exact 
trajectory  of  case  B.  The  LIN  method  gives  an  u-trajectory  that  is 
almost  identical  to  the  x-trajectory.  The  final  heading  is  about  0.05 
radians  less  than  the  desired  heading  and  the  final  point  is  about  one 
foot  away  from  the  desired  final  position.  The  best  x-trajectory  is 
obtained  by  the  SE  method  and  the  best  u-trajectory  is  obtained  by  the 
CG  and  LIN  methods. 

Table  7.4  and  Figures  7.33  - 7.37  contain  the  results  for  case  D. 

The  CG  and  the  ME  methods  give  very  good  results  for  both  the  x-trajec- 
tory and  the  u-trajectory.  The  SE  method  gives  a good  approximation  of 
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the  x-trajectory  but  due  to  the  slowness  of  building  up  a large  angular 
velocity  the  controls  take  too  long  to  reach  the  stops  and  u-trajectory 
is  not  close  to  the  x-trajectory.  The  LIN  method  gives  very  poor  results 
when  R is  held  constant  and  not  much  better  results  when  R is  exact. 

The  x-trajectory  in  Figure  7.36  exhibits  too  sharp  of  a turn  in  the 
middle  of  the  trajectory.  Therefore,  the  LIN  method  does  not  give  good 

; results  for  cases  of  long  periods  of  turning. 

i 

In  case  E the  time  interval  is  the  same  as  in  case  D but  the  amount 
of  turning  is  much  less.  The  results  are  given  in  Table  7.5  and  Figures 
i 7.38  - 7.43.  The  head  of  the  torpedo  is  the  aimpoint  except  for  the 

, LIN  method.  The  SE  method  gives  an  x-trajectory  that  is  more  like  a 

continuous  curve  than  like  the  exact  trajectory  and  its  final  position 
I is  close  to  the  desired  final  position.  The  u-trajectory  is  almost  a 

straight  line  along  the  X-axis.  The  CG  method  didn't  converge  at  all. 
When  the  initial  trajectory  was  given  to  be  a 0.4  second  curve  followed 
I by  2.6  secs,  of  straight  trajectory  the  CG  method  converged  in  one  step 

to  almost  the  same  controls.  Instead  of  no  rudder  deflection  for  the 
last  2.6  seconds  a small  deflection  was  given,  just  enough  to  almost 
reach  the  desired  final  point  with  a slightly  larger  heading.  The  ME 
method  yields  an  x-trajectory  that  turns  slowly  and  doesn't  reach  the 
I final  point  or  heading  and  an  u-trajectory  that  turns  too  much  and 

misses  the  desired  final  position  with  too  large  of  a heading.  The 
x-trajectory  and  the  u-trajectory  of  the  LIN  method  are  very  close  to 

t 

one  another.  The  final  position  is  close  to  the  desired  position  only 
the  final  heading  is  larger  due  to  the  slower  turnrate  and  hence 
longer  time  before  the  torpedo  starts  travelling  in  a straight  line. 
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The  LIN  method  gives  the  best  results  for  case  E. 

From  the  above  results  it  seems  that  one  can  formulate  the  rule  that 
for  the  fixed  time  torpedo  control  problem  the  best  results  are  obtained 
by  the  ME  method  or  the  MCG  method  if  the  torpedo  is  always  in  a maximum 
curve  but  if  the  trajectory  includes  periods  of  straight  trajectories 
then  the  LIN  method  is  the  best  one  to  use.  We  used  the  LIN  method  for 
longer  periods  of  time  but  since  these  cases  were  used  with  the  LSE 
method  for  the  time-optimal  torpedo  control  problem  the  results  will  be 
presented  later  in  this  chapter. 

Now  the  results  for  the  time-optimal  torpedo  control  problem  are 
presented.  The  only  method  which  is  directly  applicable  to  time-optimal 
problems  is  the  Epsilon  Technique.  Therefore,  the  only  two  methods  used 
in  the  time-optimal  torpedo  control  problem  are  the  SE  method  and  the 
LSE  method.  Originally,  the  geometry  of  the  problem  was  as  follows; 
the  torpedo  is  located  at  (-1000,0)  in  the  inertial  coordinate  system 
heading  straight  down  the  X-axis  and  the  target  is  located  at  the  origin 
heading  straight  up  the  Y-axis. The  aimpoint  was  the  point  of  intersection. 
For  these  runs  the  transformation  of  the  state  equations  changing  the 
desired  endpoint  from  a point  on  a known  curve  to  the  origin  as  given  in 
chapter  two  was  used.  A few  of  the  cases  are  presented  in  Table  7.6  and 
in  Figures  7.44  - 7.47.  If  both  the  initial  and  final  endpoints  are 
free  then  converges  in  very  few  steps  to  the  torpedo  heading  straight 
to  the  intersection  point.  If  the  initial  endpoint  is  fixed  then  the  SE 
method  gives  good  results  only  if  the  target  is  stationary  and  the  tor- 
pedo is  facing  the  target.  When  the  target  was  moving  the  SE  method 
generated  a good  x-trajectory  as  in  Figure  7.45  but  the  u-trajectory 


barely  leaves  the  X-axis.  Therefore,  in  an  attempt  to  understand  why  the 
SE  method  was  failing  to  generate  the  proper  controls  to  coincide  with 
the  x-trajectory  cases  A - D were  tried. 

The  results  for  case  A are  given  in  Table  7.7  and  in  Figures  7.48  - 
7.53.  Initially,  the  final  endpoint  was  fixed.  We  found  that  it  is 
necessary  to  increase  the  number  of  sample  points  in  order  to  obtain 
a good  x-trajectory  with  reasonable  controls.  When  the  final  endpoint 
is  allowed  to  be  free  the  x-trajectory  is  almost  the  same  as  when  the 
final  endpoint  is  fixed  but  it  takes  much  longer  for  the  angular  veloc- 
ity to  increase  and  the  controls  to  reach  the  proper  value.  The  ME  method 
gives  a very  close  approximation  of  the  exact  trajectory  with  its  x- 
trajectory  and  the  controls  and  u-trajectory  are  exact.  The  LSE  method 
gives  a good  approximation  and  is  slightly  better  than  the  SE  method  in 
that  the  final  heading  is  closer  to  optimal.  The  SE  method  slightly  un- 
derestimates the  optimal  time  and  the  LSE  method  slightly  overestimates 
the  optimal  time.  Comparing  these  results  with  those  from  the  fixed  time 
torpedo  guidance  problem  we  see  that  the  SE  method  and  the  ME  method 
give  the  same  results. 

In  Table  7.8  and  Figures  7.54  - 7.63  the  results  for  case  B are 
presented.  The  results  for  case  B are  similar  to  the  results  of  case  A. 
For  the  SE  method  the  fixed  final  endpoint  case  gives  better  results  than 
the  free  final  endpoint  case.  If  the  torpedo  is  initially  in  a curve 
then  the  SE  method  and  the  ME  method  give  exact  controls  and  almost  exact 
x-trajectories . For  the  case  of  free  final  endpoint  the  method  proposed 
by  Hewett,  in  reference  [10],  of  including  the  second  order  terms  in  the 
Newton-Raphson  step  was  tried.  There  was  no  improvement  in  the  results. 
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Another  modification  that  we  tried  was  letting  x^,  the  velocity  of  the 
torpedo  parallel  to  its  length,  be  constant.  For  this  case  the  system  is 
known  to  be  controllable  and  an  optimal  control  is  known  to  exist.  The 
results  are  very  close  to  the  cases  run  with  no  constraint  on  x^  and 
there  is  no  improvement  in  the  value  of  the  cost  functional  or  the  number 
of  steps  needed  to  converge.  A case  was  run  where  the  coefficients  of 
Rayleigh-Ri tz  expansion  of  the  state  vector  were  initial lized  to  their 
optimal  values  but  the  final  time  v/as  not  correct.  The  x-trajectory  is 
identical  with  the  previous  case  but  the  controls  and  the  u-trajectory 
are  much  closer  to  the  optimal  controls  and  trajectory.  The  LSE  method 
g.'ves  very  close  agreement  between  the  u-trajectory  and  the  x-trajectory 
and  both  of  these  trajectories  are  very  close  to  the  x-trajectory  given 
by  the  SE  method. 

The  results  for  case  C are  given  in  Table  7.9  and  Figures  7.64  - 
7.72.  For  the  SE  method  all  of  the  x-trajectories  are  close  to  optimal 
with  the  free  final  endpoint  case  having  a slightly  lower  heading  than 
desired  and  the  fixed  final  endpoint  case  having  a larger  heading  than 
desired  at  the  end  of  the  run.  The  only  case  where  good  controls  were 
obtained  v;as  when  the  coefficients  were  initial lized  to  the  optimal 
values.  In  one  case  the  method  of  steepest  descents  was  used  for  the 
first  nine  iterations  but  no  improvement  was  obtained.  The  LSE  method 
has  good  agreement  between  the  x-trajectory  and  the  u-trajectory.  The 
final  heading  is  slightly  larger  than  the  desired  heading  and  the 
controls  are  the  best  obtained.  The  ME  method  gives  the  same  controls  as 
for  case  B. 

The  results  for  case  D are  presented  in  Table  7.10  and  Figures  7.73 


- 7.78.  The  SE  method  gives  a good  approximation  of  the  exact  trajectory 
with  the  x-trajectory  but  the  u-trajectory  is  not  very  close  to  the 
x-trajectory  as  can  be  seen  in  Figure  7.74.  Increasing  the  number  of 
sample  points  doesn't  give  any  improvement  in  the  results.  This  is 
probably  because  the  number  of  sample  points  would  have  to  be  increased 
by  several  hundred  points  before  any  improvement  would  be  noticed.  The 
ME  method  gives  perfect  controls  but  the  final  time  was  too  short  and 
the  x-trajectory  is  not  good.  The  LSE  method  gives  a good  approximation 
of  the  controls,  the  x-trajectory  and  the  u-trajectory . 

From  the  above  discussion  we  see  that  the  LSE  method  gives  the  best 
results  for  the  general  time-optimal  torpedo  control  problem.  Therefore, 
the  LSE  method  was  applied  to  several  cases  with  longer  ranges.  The  cases 
are;  case  E with  the  initial  estimate  of  time  being  slightly  less  than 
exact  and  with  the  initial  estimate  of  time  being  exact,  the  target 
being  the  origin  with  the  torpedo  located  at  (-500,0)  with  initial 
headings  of  0.1,  0.2,  0.3,  0.4  radians  and  with  the  torpedo  located  at 
(-750,0)  and  (-1000,0)  with  an  initial  heading  of  0.1  radians.  The 
results  are  given  in  Table  7.11  and  Figures  7.79  - 7.94.  The  results  of 
the  Liri  method  applied  to  the  fixed  time  torpedo  control  problem  are 
presented  along  with  the  results  of  the  LSE  method  applied  to  the  time- 
optimal  torpedo  control  problem.  As  can  be  seen  in  Figures  7.79  - 
7.94  the  SE  method  applied  after  the  LIN  method  did  little  more  than 
increase  the  final  time  to  be  sufficient  for  the  state  vector  to  reach 
the  final  endpoint.  The  LIN  method  gives  a trajectory  that  reaches  the 
final  endpoint  but  it  is  not  the  optimal  trajectory  of  a maximum  turn 
until  the  torpedo  is  facing  the  aimpcint  and  then  a straight  trajectory 
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to  the  final  endpoint.  It  is  the  feeling  of  the  author  that  the  reason 
for  no  improvement  is  the  inability  to  use  a sufficient  number  of 
sample  points  in  the  SE  method.  The  LSE  method  is  an  improvement  over 
applying  the  SE  method  alone  in  that  the  SE  method  didn't  generate 
controls  that  guided  the  torpedo  to  the  desired  final  endpoint. 
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(a,b)  a -value  of  control  variable 

b - sanple  number  to  which  a correspond; 
Values  vary  linearly  between  points 
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TABLE  7.2:  continued 
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Figure  7.1.V  Fixed-Time  Torpedo  Control 
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Figure  7.20.  F’ixed-Timc  Torpedo  Control  Problem.  Case  B:  C(i  Method. 
Torpedo  Initially  in  a Curve. 


COMPUTATIONAL  TECHNIQUES  FOR  TIME-OPTIMAL  CONTROL  OF  TORPEDO  DY--ETC(U) 
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TABLE  7.3:  Case  C Results 
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TABLE  7.3:  continued 


Figure  7.24.  Fixed-Time  Torpedo  Control  Problem.  Case  C;  CG  Method. 


TABLE  7.4;  Case  D Results 


TABLE  7.5:  Case  E Results 


Figure  7.39.  Fixed-Time  Torpedo  Control  Problem.  Case  F:  SF  Method,  Aimpoint  Head 
of  Torpedo,  x-Trajectory. 


Fixed-Time  Torpedo  Control  Problem.  Case  F:  CG  Method,  Aimpoint  Head 
of  Torpedo,  Initial  Trajectory,  0.4  sec  Curve,  2.6  sec  Straight. 
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TABLE  7.8:  Case  B:  Tine  Optinal  Torj,'edo  Control  I^blem 
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TABLE  7.8;  continued 
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TABLE  7,9:  Case  C:  Tiine  Optinal  Torpedo  Control  Problem 


TABLE  7.10:  Case  D:  Time  Optimal  Torpedo  Control  Problem 
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Figure  7.76.  Time-Optimal  Torpedo  Control  Problem.  Case  D:  MF  Method,  12  Basis  Functions, 
50  Sample  Points,  Free  Final  Fndpoint. 
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TABLE  7.11:  LSE  Method  Applied  to  Long  Ranges 
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TABLE  7.11:  continued 


TABLE  7.11:  continued 


7.11:  .jntinued 


Figure  7X1.  Fixed-Time  Torpedo  Control  Problem.  Case  F UN  Method,  A = 60.000. 
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Figure  7.83.  Fi.xed-Time  Torpedo  Control  Problem.  Initial  Range,  .SOO  ft.  Initial  Heading 
0.1  rad.  LIN  Method,  X = 60,000. 
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Figure  7.84.  Time-Optimal  Torpedo  Control  Problem.  Initial  Range,  50  ft,  Initial  Heading, 

0.1  rad,  LSF.  Method,  iJ!  Basis  Functions,  51  Sample  Points,  Free  Final  Endpoint. 
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Figure  7.8S.  Time-Optimal  Torpedo  Control  Problem.  Initial  Range,  500  ft.  Initial  Heading, 

0.3  rad,  LSF  .Method,  12  Basis  Functions,  51  .Sample  Points,  Free  Final  Endpoint. 
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Figure  7»0.  Time-Optimal  Torpedo  Control  Problem.  Initial  Range.  500  ft.  Initial  Heading, 
0.4  rad,  LSI;  Method,  12  Basis  Functions,  51  Sample  Points,  Free  Final  tndpoint. 


160 


1 


M 


Chapter  Eight:  Conclusions 

In  this  dissertation,  several  computational  techniques  are 
applied  to  the  torpedo  control  problem.  The  SE  method  gives  a good 
x-trajectory  but  the  u-trajectory  is  unsatisfactory  in  many  cases.  The 
ME  method  and  the  CG  method  yield  good  results  for  cases  A,B  and  D but 
not  for  cases  C and  E.  The  LIN  method  yields  poor  results  for  cases  A, 
B and  D and  fair  results  for  cases  C and  E.  The  LSE  method  gives  good 
results  for  all  of  the  cases.  There  is  close  agreement  between  the  x- 
trajectory  and  the  u-trajectory  and  the  final  endpoint  is  always 
reached.  Therefore,  the  feasibility  of  using  the  LSE  method  to  solve 
the  torpedo  control  problem  has  been  shown.  It  is  possible  that  the 
LIN  method  is  not  the  only  method  that  could  be  used  to  initialize  the 
trajectory  of  the  state  vector  in  the  Epsilon  Technique  to  obtain 
adequate  results. 

The  computer  limitations  of  limited  computer  storage  that  was 
discussed  at  the  end  of  chapter  three  severely  limit  the  optimality  of 
the  results  of  the  SE  method  and  the  LSE  method  to  the  time-optimal 
torpedo  control  problem  for  the  cases  where  the  trajectory  is  of  a 
time  duration  greater  than  two  seconds.  It  is  our  conclusion  that  the 
results  would  be  closer  to  optimal  if  a sufficient  number  of  sample 
points  and  basis  functions  could  be  chosen. 
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Appendix  A:  Computer  Programs 
In  this  appendix  a brief  description  of  the  different  main 
programs  and  subroutines  is  given  followed  by  a listing  of  the 
computer  programs.  In  most  cases  the  names  of  the  variables  that  were 
used  in  the  text  are  used  in  the  computer  programs. 

MAIN 

This  is  the  main  program  used  for  applying  the  Epsilon  Technique. 
The  order  of  the  program  is  to  define  the  initial  conditions.  Then, 
the  epsilon  cost  functional,  its  gradient  and  Hessian  is  computed. 

Then,  system  (3.13)  is  set  up  and  solved  to  determine  the  new  vector 
A.  Finally,  the  epsilon  cost  functional  is  again  evaluated  to  determine 
if  the  method  has  converged.  If  not,  the  iteration  is  started  again. 
Subroutine  FCT 

This  subroutine  is  used  by  MAIN  to  calculate  the  state  vector, 
its  derivative  with  respect  to  time,  and  the  derivatives  with  respect 
to  A using  the  Rayleigh-Ritz  expansion. 

Subroutine  DFCT 

This  subroutine  is  used  by  MAIN  to  calculate  f(t,)^,u)  and  the 
partial  derivatives  of  f(t,^,u)  with  respect  to  x.  The  optimal  value 
of  u is  calculated  in  this  subroutine. 

Subroutine  COST 

This  subroutine  is  used  by  MAIN  to  calculate  the  cost  functional 
and  its  derivative  with  respect  to  the  time  step. 

Subroutine  RDLU 

This  subroutine  is  used  to  calculate  the  solution  of  the  equation 
Cx.  = by  a LU  decomposition  and  back  substitution  using  double 
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precision  arithmetic.  This  subroutine  can  also  be  used  to  compute  j 

the  inverse  of  C.  This  subroutine  is  used  by  all  of  the  different  | 

methods.  i 

CONGRA 

This  is  the  main  program  used  in  the  conjugate  gradient  method. 

Steps  I - 6 given  in  chapter  4 are  implemented. 

Subroutine  FCT 

This  subroutine  is  used  by  CONGRA  to  calculate  f(t,^,u)  and  its 
derivatives  with  respect  to  x and  u. 

Subroutine  COFCT 

This  subroutine  is  used  by  CONGRA  to  calculate  all  of  the 
variables  given  in  Step  3 except  for  and 
MNLNRL 

This  is  the  main  program  used  to  apply  the  method  of  lineariza- 
tion about  a known  trajectory  as  given  in  chapter  5. 

Subroutine  COEFC 

This  subroutine  is  used  by  MNLNRL  to  calculate  the  constants 
used  in  the  equations  of  motion  of  a torpedo. 

Subroutine  FC 

This  subroutine  is  used  by  MNLNRL  to  calculate  f(t,^,u). 

Subroutine  FCEV 

This  subroutine  is  used  by  MNLNRL  to  calculate  the  Riccati 
equation  for  the  method  of  linearization  about  a known  trajectory. 

MNLNEP 

This  is  MAIN  modified  to  call  MNLNRL  first  and  to  use  the  output 
of  MNLNRL  to  initialize  the  state  trajectory  in  MAIN. 
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00  6A  X ■ 1 , 9 

0i»040(LiJtX)..DPrX|J.x*3)*0xD40ix»2I.S3C'E 
46  lPIJ,E0.x»2IO>.O4oiL,J,xl«CwO4c(L,j,X)xSIOE.DXOC*a(JI 
COnTinlE 

00  8 7 J"  1 ,N 
OPOT 1 J I *0. 

00  87  1 • 1 , S 

87  :pOT(J1p0FDT(ji.oP0xij,I  i.DxOTl  i ) 

OC  J J«  1 . b 

«lL.j)*<XO(JI-PIJIi«S30f 

3 D«007Il.J'««ILiJ1/I2.»0^1»IOXOOT(jI-oP07|JI  IaSOOE 
IP|L,LE.NpSITI“E.IL-I  , l•L)7/^PS 
IP(L,G7.NPS|9I“E.lL->.PS'»rT 

PRINT  10b,TI-E,IXIIl,|.l,6l,IXOlII,I.l,il,|P(n,I.|.Al 
10b  P0R“4T(4,»  T . ,PlO,6,9>->  • ,6E15.5./,6m  xO  • ,6E2C.5i/iSN  P ■ , 

% 4E  2;.s I 

XT»Bj.xT4CG»x70»0T 
r T A P G • 7 7 4 0 G • * 7 0 • 0 T 
I C 0 N T ; U E 

PRINT  50  , (I  ,CEL I I I , I . I . Np  I 
50  P0R"aT(5h  I . ,I3i2x,6hOel  • ,E20.1bi 

call  COST  I ?T , ^P , j ,cG0T , X , XPP , n I 

00  9 J. 1 , ' 

00  9 |.l 

9 A 1 ( J.9-9.  n .A  ( J , I I 
00  67  1*1,9 

41(N.“.I|.4C0II»2I 

6’  continue 
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NOSC-TR-124 


V 

I 


11 < 102  I ‘OT 
COEF.COEFO 
DO  1 *t»  I • 1 , I 02 
1«»  M0(  I I**!  ( I I 
00  S J«1 
DO  5 L»1 

5 Ml  ( J«NF.NF»LI«m(1.iJ) 

Ml ( 101 )■& 

DO  * J«l 
00  4 t»l 
00  7 F» 1 , N 
00  7 I»1 

7 DM1D»1(J»NP.nP*L,«»m.m*1|.D»D1(L,Ji*,II 

00  *B  1 , H 

48  DMlDll IJ»NP.NP*L,S»N»F|«0W010(1,,J|FI 
4 DM1 D* I IJ«NP.NP»L , I 02  IPOpDOTil.J) 

DM  ton  I lOI  |I02  I •OGDT 
DO  9 I • 1 , I 02 
DO  8 J«1 , 102 
HQRh I I , J I aQ I 
DO  8 Ka I , 1 D 1 

8 MGRM|l,J|aHGRMlI,J)*DMiO*iu,ii,DHDil(K;j) 
MOPC I 1) ao . 

DO  9 Fa  I , to  1 

9 HDPC(  1 laMDHCl  I laOM|D*l iF, I laM)  (F  ) 

CVaO. 

DO  12  t*l • lOl 
|2  CVaCv*Ml  I I I*a2 
CVOaCV 

PRINT  I 02 , ILOOP ,CVO 
102  F0Rn*T(5h  I a , IS i2X ,E2°. 1 5 I 
DO  !•! ,102 

3<4  HGRM  ( I I I ) aHGRM  I I , 1 I *D  I iG 
DO  902  lal , ID2 
DO  902  ja I , 102 
902  HGRMI ( I iJIaMSRHl I , J) 

CALL  TREOI  I 102, ID2,mgR„I ,D,E,E2I 

CALL  I nTOL 1 I I D2 ,0 ,E , IER" I 
PRINT  900 , lERR , ( 0 M ) , 1 a 1 , 102 ) 

900  FORnaTI 7MIER  a ,IS,9HD  • , 1 00 ( T i e , E 2o , 1 5 , / 1 1 

call  RDLU t ngRm , mgRm I , md"C I HOnO , I 02  ,0,01 

GO  To  74 

75  DO  79  I«1 ,102 
7"  MDND(  I laMONCt  I I 
74  continue 

DO  10  ■■! , 102 
A I I I I aA  H II aCOEFaHDMD I 1 > 

10  continue 

00  II  J9 1 ,N 
00  1 I I ■ I ,N 

11  A I J,  I I aA 1 ( JaM.a* I | 

print  lOJ, ( ( A( I , Jl , J,1 ,P| , i.i ,ni 
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I 


1C3  rO»HAT(*(lH  ,*E20.15,/,f3,AEr0ti5./)',//l 
00  69  !•!  ,*t 
*00  I I »2 I •* 1 I » I I 

continue 

DTa* I 1102) 
inOT.LT.O.IOTa.l 
1M0T.ST,.5)0T..5 
CV.O, 

SOCE.SQST I OT/E^S ) 

00  151  L»1|NP 

C*LL  ECT(X,*D,»,*0, *00, 0x0*, 0*00*, 0«0*0t 0*00*0, OXOT,DXOOT,l,N>,NP 
f ,0T , *L^H* , *or ) 

C*LL  DrCT(F,*,*0,0r0X,L»N,0T,0tL',N»’,*T0,YT01 
DO  151  J>1,N 

Ml L I J I ■ I *OIJ I >F ( j , I *500^ 

CVmCv*m(L, Jl*«2 
151  continue 

C*LL  COST ( OT ,NP ,0 ,OSOT ,* , XfP ,N ) 

CV«Cv*6»6 

IFKv,LT,CVO)GO  TO  153 
COEF.COEF/H. 

DO  15S  I>1 , 102 
ise  *t ( I i«*io(  1 1 
60  To  76 
153  continue 

print  IOM , I *0 I I I , I • 1 ,N ) 
print  1 Or , I *00 1 I ) , I ■ 1 ,N I 

10»  FORH*T| IH  ,3E20. |5l 
DE>Cv6*6 
TF.InP-1 1»0T 

1 LOOP* I LOOP* I 

print  1 00 , I loop , TF ,cv ,CE 

100  F0Rn*T(5M  I ■ ,Ir,2*,5nTf  , , E 1 5 . 5 , 2 * , 5NC V • , E 1 S . 5 , 2 * , 5MDE  • i 
lElS.SI 

DO  UO  1*1,102 
1*0  *15  11, ILOOP )■*  1 I I 1 

IFl Il00P.lt. 30  160  TO  2° 

00  1*2  1*1 ,2S 
DO  1*1  J* 1 , 102 

1*1  *lRIJ,I)**lSlJ,II*M*|Slj,i«||.*|Sij,l|la«2i/|*|S,j,|*2|. 
f 2,**1SI J,  1.1  |**iS( j, I ) I 

DO  1*3  *• 1 ,N 
00  1*3  L*1  ,M 

1*3  * I X ,L  •**  IRt  **►'•". L . I ) 

1*2  PRINT  103,1  (*!*, LI  ,L*I  ,'■1  ,«*1  ,N| 
print  145, I * iR1  102, I ), 1*1  ,281 

1*5  F0Rn*T17m  *1P  • , 1 T8 ,E iS, 1 0 I ) 

TPP*,01 

T*0. 

DT*Tf/I 3T7, 

00  120  1*1, N 
120  XII  I **0 I I I 
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I 


L«NP.l 

00  121  l»l  I 1377 
I r ( 1 ,CT . I TEST IL«L» I 
inL,LE«NP*NPSI  ITEStO.3 
IFiL.6T.NP»NPS)ITESTD»iSj 

IF(I,GT.ITESTMTCST,ITeSt»1TCS70 

CALL  OrCT(F,x,XO,Orox,i.»SiOT,OtL,NF,xTO,rTOI 

00  122  J* 1 |N 

122  XI JlaXl JI*0T<F( J| 

T.T*0T 

IFlT.LT.TFPlGO  TO  121 
TPP»TPP*.01 

P«IHT  123,T, (*( j| ,j.i ,nl 

123  FOHP»T(SM  T . ,FiO,3,*E>5. 101 

121  continue 

call  OPNPLT 
call  XTEPLT 

call  linpltixplt,ypi.t,n^i 
call  curve t XPTO , yPTG.NP .0  I 
call  ENOPLII) 

GO  To  21 
230  continue 
STOP 
END 
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I 


FCT 


sue^ouTiNf  fcTi«,«c,A,*c,4ooiP»n»,o»oo»,oxo*o,D»oo»r.c»oT,DiODT,L, 
«N,M,»JP,DT, ALPHA, aOF) 

DIHEnSIOS  «IN),ApiN),Dxf>AiN,H|,|3»COA(N,H),OXDAO(N),rXOOAO(NI, 

*DX0T(N) .OXDOTIM ,A|s,N) iAOIS) ,Ano(N» 

» , AOF I N ) 

oefimE  c ( I ) .s  I “>  1 1 »p  I • t /Tf  ) 
define  00(  I I*I»PI/TF«Cc5(  I«PUT/tFl 
PI«3. 1N15926E36 
BL»L 

T.PT.(l-I  . I 
TF»InP«1.)«DT 
DO  I J»  I , N 
X ( J I aO* 

XC ( J ) >c • 

DO  2 lal.H 
0X0* I J,  1 I «0 I I I 
DXDDAlJ.MaCOlII 
IIJIaXIJ)*A|j,I).DxDA(j,|l 
2 XO ( J I »x D I J ) ♦ A Ij , I I t 0 XDC* IJ , M 
X(JlaAO(JI*AOD(JI*T/TF.»|JI 
XO ( J I aAOO ( J I /Tf  axD I J ) 

OxD*cIJI«T/^F 
OXDDaO(JI"1./TF 
OXDT ( J ) aO . 

OXDDT ( J I a.xD IJ  f /OT 

1 continue 

bETUBn 

ENC 
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I 


c 


SUBSOUTiNf  DrCTIF,r,irO,DfO»,L,f-,DT,OEL,‘'P,»TO,'rTD| 

OIMEnSIOM  riN),X(K),iiD(N),DEOx(N,NI,OCLl»<P| 

COH«O^J/*/I'FDn»  16,6,6  I iPUTEST 
RK»2S« 

X la  • 9 
RMOaJ, 

RLa I 5« 

0 B * X a 1 6 . / 5 7 . 2 9 5 7 8 
XTla.7,5 

RXa 1 ,C7/Z> 1 
DP  a • J 5 

RFaPL/l  I • IZBaSRPTI  11)  ) 

RXla,65i«RF«a(»l  tH9) 

RKZa , 73«RF . 1 1 39 
RHLasOO  . •(  laPKl  )/37,2 
RMTa800.PI l♦Rx2)/37,2 
»XPa,‘*2*RF**.331 
RJZaMOO.P3HO.PPKP 
RLTlPaZ.S 

THRUSTa(RNO/2.)PlXPOPp6^.56P67,56 
M8a)RH0/2.  IpIIpRlTIP 
eiaC0SIX(6l I 
B2a>S  I N I X ( 6 I ) 

B3a-02 

BMaP  1 

RMOMaRHO/2 . 

H 1 aPHOHPllPDP 
HZaPwOMPll 
MOaPwOH. 1 laSL  P . 3 
M6a>) .3p|RMC/2. IPllaPL 
M7aPM0M»llaRLPPl-Pl*.2) 

I 0Ca300 

1 OCaiOO 
Ka). 

I P ( L , GT  , NP  ) K aL  -►iP 
IFIL.G7.NP) GO  70  5 
IF ( L.tO.O IGO  70  2 

laxD)MI»x(3lp|XI5lp|M3.RML)-M7.)()a)-  Mg,|xlM)»XTlPl(5l ) l/R"7 
0a>  m8p»K p X I 3 I a X ( 3 ) /B“T 

CaXD(5,t-(M6PX(3)pX|M)»H7.X(3laXI5)«  M8.xTlPX(3lptX|MI»XTlaxlSI)l 
* /P  JZ 

Daa  M8pRKPX71pX  I 3 ) ai I 3 ) /R jZ 
OCLU  Ip'I  iPBaCPOl/IBPBaOaDI 

IMlBSlDEL(X)).GT,DUTES7)0EUXIaS16N)0Mlx,DtL(X)l 
5 C0N7lNUf 

v7aX(«)P171alI5).RK,x(j),[)EL(|,) 

GO  70  3 

2  IFlNP,Lr.lOCIV7a*|a).xT*Pll5)»Rr.X(3)PD“lX 
IFtNP,G7.I0ClV7al)*),xT*PX(5) 


3 continue 

FIiI.B;»«I3I»P2»*Ih|.«tO 

riZI«P3*«(3l»P‘*»«l‘<l-TTO 

rl3l«IP"T.X|Ml««(S)»TMROST-Ml»XI3l»F(3|l/RKt 
riM)»XI3l»IIIpl*|H3-RMLl“H2««l*.|-  h8»VT)/R"T 
F(5l.|M6»»l3l»«l*'l»HT.il3l.»(5l.  XT*»"e»Xl3l»VTl/Rj2 

r I 6 I ■ X I s I 
IFIl.GT.npibeT’jpn 
IFIL.ES.O  IBETuBn 
CFCX I 1 I 3 I >8 1 
DFDX ( 1 , N I ,B2 

OFOXI  1 ,6)»B2*xl3l-Bi»i|'*) 

DEC  X I 2 I 3 I *8  3 
OFDXi2,h|.Bi 

CFCX  ( 2,*>«ei«xl3l»B2»Xl‘'l 
OEDX|3,3I»-2.*“1*XI3i/B"l 

DECX|3,••I•XIBI•P"T/<)N^^ 

DFOX|3,PI»xi‘*I*P'‘T/rnl 

OFCXi*i,3I»FihI/X(3)*  Ng,BK,pE|_(Ll*xl3t/P>*T 
CFDX|M,>*).|-n2*XI31-  Ng.,|3))yPKT 
CFD»I‘*iPI*IX!3I»I’‘3-Pkl’*  M8»x(3I»»Tai/IJmT 

CrcxiS,3l«(N6»X(<<).«7.x(&l»  XT»»H8»(yT.X(3l«»X»0EL(LI I l/Pjr 

Crcx(5,'*l»IH6»x(3).  X^**h8»x(3II<''’j2 

0rDX|5,5l«(M7»xi3l-  x^‘»nb»xI3I»xTi|/Oj2 
CFDI 1 6 I S 1 ■ 1 • 

OEDCx ( I . 3 , i I .82 
orccxi 1 iM,6i.-Pi 
CFDCxI 1 , 3 ).?2 

crcox I 1 .6,8 1 .-B  1 
CFOOx I 2 , 3 , 6 I •?  I 
OFCDx 12,8,61.82 
CFDCx ( 2 , 6 , 3 1 .P 1 
0FC0XI2, 6,81.82 
0FCDx(3,3,3I.-2.»81/9ml 
orccx  I 3 , 8 ,s  I .b“T/rmi. 

CFDCx I 3 ,5 , 8 I .R"T /Bml 

3FD0xl8,3,3l.8p»Bx.DELI''l/PPT 

CFCCxl8,3,8).(-H2-Ng|/BMl 

0rCDXI8,3,'il.|H3-BMl.-M8*XT61/P8T 

0FDDx(8,8,3I.I*«2*M8)/b»'T 

DFC0xl8,s,3).(,'3.BML.«g.xT*)/p«T 

0ECDxIB,3,3|.2.*xT6.h8»Px»DELIX)/BJZ 

CrcCxlP,3,8|.l«6-XT6.M8lxBJ2 

CECOll5,3,fl.lN7.XT».xT6.M8)/RJZ 

OFDDxl5,8,3l.l«6-XT6»Me'/BjZ 

OFODxI‘,,S,3I.|h7-xT4.xT*»h8I/bj2 

RETURN 

END 
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I 


I 


5ue<»0UT!Nf  C0Sfl0T,KP,Q,D60T,I,xrP,N) 

0 InESS I ON  X I N I , IFP I N I 

G«SQPT ( I NP- I . I POT  I 

DGDT. ( NP- I . I / I 2. pGI 

RETURN 

END 


pDlu 


SUB"OUTInC  R0LU(»,*|Nv,CvECTB,X,NC*'’,M0fA,|NVr6) 

OIHEnSION  •lnC*P,NC*P|,*lNV(NC*P,NCiPl,CvCCTR(MC»Pl,*|NC*^) 
double  PBECISITN  *,*INviCVECT9,x 
double  precision  OETR .oCono ,0ET I 

double  precision  H*xPIvtFLTTS,TSC,P|NlNUSI,NTSI&N,OfHKH*,OCMKPN 
C • • • • • pROGRXhmED  By  don  6*R*Cm  nuC-SO  code  601  • • • • • 

C . . . • M1»PIV,"INUSI .nEsIGN  must  be  typed  REIL  FOR  SINGLE  PRECISION 

DIMENSION  ROW ( 20a ) 
integer  ROM 

COmmdn/v*lOET/OETR,OET1 
COmmoN/DECOND/DCDNO 
OET I«OiO 


solution  of  TmE  STSTEm  tX.B  for  I (InvFGrOI 

c OR  compute  * inverse  (INVFGpII 


M I NUS 1 •* I . 0 
NTSIGN*!  .0 
MDFsMQrA 
L 1 M0«  i 

I F ( I nvFG . E3. 1 I LIMD«NC*P 
DO  S2S  I !• 1 , L Imd 

ifiimvfg.es.oi  so  to  22S 

DO  MIS  <K» I , NC  »P 
MlS  CVECTRiFKIrO.O 
CVECTRI  I I ).l  .0 
mDFr  1 

IF  n 1 .ES,  I V MOF  mo 
225  IF(MdF.EO. II  GO  TO  fcS 

C • • • • > DECOMPOSE  THE  matRi*  » INTO  LOmER  AND  UPPER  TRIANGULAR  MATRICES 
c pMERE  AmLU  (miTM  Column  pivcTInGI 

C • • • • a the  original  a matrix  IS  destroyed 

DO  1 lali  NCAP 
I ROM  I n a I 

MAxP I Va. 1 ,0 

NP I V OT  a 1 

00  M I a I , NCAP 

FLTTsaOABS I A I I , M I 

I F I MA XP I V . GT , FL T»S  I GO  Tq  m 

MA  xP I VaFLTTS 

NP I YQTa I 

M continue 

PRINT  7TR 
77»  FORMaTI IH| I 
L ■ i 

7*1  FORMAT! 13, 2X,*mmaxpiv  ,,015.61 
IF  1 MAXPIV  ,EG.  0,  I 6°  TO  625 
DCmxmXrMAXP I V 
OCmkmNpMAxP I V 
I F ( NP I vOT ,E0. I I GO  TO  7 
NTS  I GNaNTS I GNaM I NUS I 
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I 


ME"P«*Owl  1 ) 

ROW ( J I »POW ( NP I VOT  I 

ROW  I sP I VOT  ) «NTt  MB 
J.NP I VOT 

00  6 I « 1 , NC»P 
TSCbI II  • I I 
*I1.||bR|JiII 
4 * IJ , I 1 bTSC 
7 p« 1 ,000/* 11,11 
DO  6 ImZ,  MC*P 

e *11,11  bp.* 1 1,11 
00  55  I.B2,  *<C*P 
m*»PIVb-I,0 

ic.L-1 

00  25  1*1.,  MC»P  j. 

00  1 5 J* 1 , * 

*lt,LlB*II,L)-*IJ,LIB*|I,JI 
rLTTs.0*B5 1 * ( I ,U  I 
1 r ( Mi «P I V . GT , FLTTS  I GO  Tq  25 
iiPIV.ri.TTS 
P 1 V 0 T • I 
25  continue 

IF  I M*VP!V  .EC.  0.  I gI^  to  625 
IF(MiXPIV.GT.OCHrMx)  DC^xMx.Mixpiy 
lF(MiXPIV.LT,OCHKMNI  OC^XMN.MIXBIv 

I F 1 NB I VOT . EG . L I GO  TO  2®  I 

NTS  I GN.'-TS  I GN*M  I NUS  1 j 

NTEMP.POw ( L I 

ROW  I L I B®OW ( NP I vCT  I j 

ROW  I nP I VOT 1 .nTEmp 
J.NP 1 VOT 

00  27  LI,  NC*P  I 

T5C*  I L , I I 
*IL,II.*IJ,II 
27  * (J,  ! 1 .TSC 
2#  P.1.0  /*IL,LI 
X .L»  1 
“.L*  1 

IFIL,EG.NC*P|  go  To  55 
CO  *5  LX,  NCiP 
00  35  J. 1 , M 

35  * I L , I L*  I L , I I -»  I L , J L*  ( J , 1 I 
.5  1{L,II.P»*(L,II 
55  continue 

OCONo.OLOGIOiOCHxmxi.OlOsIOIOCMxmnI 
DETR.* I 1,11 
00  205  L2,  NC»P 
205  OtTR.DETR.*!  I , I I 
OETR.OETR.NTSIGN 


C • • B » • USES  * right  h*NC  side  B *ND  SOLVES  FOR  THE  UNXNOwN  x BT  TWO 
C SUCCESSIVE  BiCX  SUBSTITUTIONS  lT.B  *ND  UX.T 


*5  DO  *6  LI,  NC*P 
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66  I I I I aOtO 
JaROi.1  1 ) 

«(1).CVECT(»(J|/*I1,1I 
DO  85  1*2|  NC*P 
J«ROt.{  I I 
IC«I-1 

DO  ?5  L»  I , K 

7S  X(IlaXlt)*A(I,LI*XILI 
85  X ( I I . I CVECTR IJ I -X ( I I I /* ( 1 , 1 I 
XaNC»»- 1 
DO  ID5  I»1  , X 
J.NC  »p- I 
x.NC  »P 

DO  95  L*  1 , I 

XIJI»X(JI-XHKI»*(J,M| 

95  M.M« I 
105  continue 

in  lNVF5.Ea.CI  50  TP  5j5 
C • • • • . LOADS  T«E  INVERSE  OF  the  mATRIx 
DO  920  LL* 1 , NCAP 
920  AlNViLLiMlailLLI 
525  continue 
RETURN 

625  print  626,  L 

626  FORPlTl*  SINGULAR  "xTRiX  . . COLU****'  , 19"  • 6LL  2ER05*  1 

STOP 

END 
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CONCK* 


f«R<MCTfR  N«6|MPa30| 

OIhEmSION  DrOX|N,NiNP),*(N,NPI,OfOUlN,NP|,DEL(NPIiP(N,NP|, 

*PnN,NPI,G(KP),OHX«(N,Nl,DHxU(N,NP),t((NP|,r{N,NP),C(N,NPI,0(NP|, 
*cr I N,NP ) ,5  I NP I ,DEL I ( NP I , *lP| 3 I .OEL?( NP) ,81 J.3 I .COEM  3 I , VEL I 3 I 
CIHEnSION  VV*U 10) ,81 I3i3I 
DIMENSION  VXU  20  I ,DXU  20  ) 

DIMENSION  F ( N,NP I ,7F ( N,Mp 1 , vr ( NP I ,V(NP ) , VSF I NP I ,SSF I NP ) , VS  1 I NP ) , 

• SSI  I np I 
NP 1 »NP-  I 
TF>,« 

TEp  1. 

OT.Tr/ ( NP. 1 . I 
E I • I 00, 

E2- I 00, 

0MiX.l5,/52,2?S78 

X I I ,  1 1 •' 1 000, 

X I 2 , 1 I pO. 

X I 3 « I I >8  7 . S6 
X I •( , I I pO  , 

X I S , 1 I pO, 

X 1 1 , I I >0, 

DO  32  1 • I ,S0 
12  OElUIppOMXX 
DO  3l  !•! ,NPl 
LPl»l 

CXLL  FCT(r,DFDX,X,N,I,pFpU,NP,OEU 

DO  2»  JP| ,N 

X I J,L )pX I J, 1 ) »OTpF Ij, I I 

2*  continue 

80  TO  31 

50  DO  S|  JP| ,N 

51  X(J,L>M»IJ,ll»0Tpl65,PF<J,n*SV,PF(J,I.l|*17,pF(J,l.2l 

• -V.pF  I J,  1-3)  )/2‘*, 

CALL  FCT (F ,OFDX,X ,N,L,dFdu,NP,DEL  ) 

00  30  J»l,N 

XIJ, L)»*(J,II*DT.(».»F(J,UPl9,PFtJ,I|pS.PFtJ,l-II»F|J',  I-2ll/2‘), 

30  continue 

31  continue 
DNp8, 

X6Fpx(*,NP I 
X I FpX I 1 ,NP 1 
X2FpxI2,NP| 

XIFpxI 1 ,NP)»DNpCOS|X*F) 

X2FpxI2,NP)pDNpSIN(x*F| 

X3FPXI 3,NP) 

XNFpX(‘),NP| 

X5Fpx(5,NP| 

print  I00,(0EL(I),X(|,iI,X(2,II,x(3,1I,X(*(,H,x(5,I),XU,1),I, 
t|pl ,np,Si 
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I 


TTI.,6 

Tn»j. 

Trp.-.oi 

DO  52  M. I , 2 
TF«TFI*l"*lil»TrO 
DO  33  I**t0,50 

33  0EL< I I»0. 

T0«0, 

00  IS  I • I ,NF 
I 8 H I I ) • I • 

1 t.00p»0 
DT.Tf/ ( up-  1 . I 
PPINT  106,M,TF|DT 

104  FORHiT I IHl  ,5M  P ■ ,IS,S*,SMTF  ■ , t 1 5 . I 0 i 5 * , SHD T • ,E16,I0I 

00  Sx  !•! ,NP 
54  5 I I I «0t 

bet»o« 

1 5 DO  It  I • 1 , N 
I 4 2 ( I • 1 I *0t 
V41.0a0« 

DO  I I»l ,NPl 
1 I>1*1 

CALL  FCT(F,DFO*,*,N,I,t)FDi,,MP,OEL) 

DO  3i|  J>  I .N 

ZF  ( j , I KOFOUIJ,  I I •«(  I I .5  ( 1 I 
DO  3<t  F»1  ,N 

34  ZF(J,II»ZF(j,Il»DFO«(J,«,lt«Z(F,n 
DO  3a  J« I |N 

I ♦OT»F  ( J , I I 

34  Z (J . I II •: I Ji I I -DT-ZF (J, I I 
GO  To  I 

35  DO  37  J»1  ,N 

*(J,in»X(J,n»0T*l55..T|J,n.5».#Flj,|.i)*37,.F(j,I.2l-P,.F(J,l-3 
SI I/Zp. 

3T  ZUiH  l-ZIJ,!  I«DT»(55.»ZFIJ,I  l-5Pt»ZF(J,I.l  l*37.»ZFIJ,I.Z). 

S». »ZF ( J, I -3  I I /Z4, 

call  FCT(F,DFDX,K,H,II,Ofdu,NP,dELI 
DO  38  J«l  .N 

ZF|J,IIlPDFDUIJtIII*w||I|«SlIII 
DO  3g  F«l  ,N 

3*  ZF|J,III»ZF(J,ni»DFD«(J,IC,in»z(Ftln 

DO  39  JP I tN 

«(J,IIl9*(JiM»0T9(9,,plj,J||»i9,,F(J,II.5,9F|j,I.ll*Ftj,l-ZII/Z4, 
34  ZIJil!!9Z(J,n*DT9|9..zE|J,ni.i9,9ZFIJ,II.5.9ZF|j,|-|l.ZF(J,I-2l 

1I/Z4, 

1 CONTiKUt 

PRIKT  100,IDEL(ll,Xll,ll,*(Z,M,XI3iIli*(4,II,l<(5,n,*(*,n,l, 

SI»1 ,nP,5I 

VAL0.El»(X(l,NPI.XlFI»»Z*tZ9(X(Z,NPI»XZFi«,2 

VALOpC|*IXII,NFI*ON9COS<XI4,NPI|.X1F|99Z*C29|XI2,nPI9DN9$|NIXI8.NP 
tl I-XZF199Z 
VAL0aVAL0*TF 
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yVALlHIaVALO 

C( 1 ,NP)«2.»Ei»ZI 1 ,NP| 

C(?,NPI»2.»t7«2l2iNP| 

Cl  3 »nP I “O* 

C M iNP I *0 . 

Cl S I NP • 

Cl Ainpi»o. 

PI|iNPI*2>*E|*IAII|NP|.A|P| 

Pll,NPI«2t»El«IXIl,NP|»*ir*DN»C0SI*l4|NP)ll 

PI2iNPI»2.«E2»I«I2,NP(.*2F| 

PI2,NPI*2.*E2*IXI2,NP|>X2E*0N<S|N|XI4,NP|I| 

PI Jinp>»o. 

PI<t,NPI«0. 

PI5iVPI»0» 

P I A I NP I «0 . 

PI4tNP)«0N»IPI2,NP|»C0S*XIA,NPn»Pll,NPI»S|NIX(4,NPnl 

DElInPI-OElINPI  I 

caul  rCTir,DrDX,X,N,NP,OF0U,NP,0tLl 
ITE-D 

VSI inpi*o. 

SSI InPI>0> 

00  » I-NPI ,1,-1 
L»l*l 

CAUL  COrcT(X,P,Pr,6,OMX*,DMxU,2,C,0,cr,N,L,NP,DrDXiOrOU,W,S.OELI 
VSFIlI**<  n*CIL)*GIL) 
ssr I L >■"1  M «S I L >*0| L I 

DO  3 J>l  ,N 

PlJi  I ) •'  I Oil.)»OT»Pr  U,L  ' 

3 CIJ. I t»Cl J,L  »-Ot.Cr ( J,L I 

ysi 1 1 iivsi ili>dt«vsfili 

SSI  I 1 >»5SI  I L t»01  »SSF  I L ) 

60  To  A 

40  00  4l  J*1 ,N 

PIJiII»P|J,Ll*0T.I65,,pF<j,ul.54,«PF(J,U*l|*37»PFiJ,U»2)- 
I4t»PFI JfL»3l 1/24, 

41  CIJ,II»C(J,L)-DT»I55,«cF(J,l)-5»,»CF(J,L»1|*37,»CFIJ,U*2I. 
S«.*CfIJiL*31  1/24. 

VSIIII*VSIILI-0T.I55.»v5FILI-54.»VSF(L»ll*3Tt«VSF(L»2l»».*VSF(L*3l 

•1/24, 

SSII:i»SSIIU-0T»l56.»5SFILI-54.»SSF(L*l|*37,»SSF(L*2l*9.*SSFa»JI 

•1/24, 

CAUL  C0FCT(X,P,PF,G,0MX*,DHXU,2,C,0,CF,N,I ,NP,OFOx,OFOU,XI,S,DELI 

SSFI I I-Xl  I-l  MSI  I I.DI I I 

VSFI  I ).WI  I-lMGI  I MDI  II 

00  4j  J»1,N 

PIJ.IM'’(J,LI-0T.I4,.PFIJ,I)*19,.PF(J,LI.6,»^F(J,L»1 )*9F(J,L*2) I 

»/24. 

42  CIJ,lMCIJ,LI-0T»l»,.CF<J,II»lf,»CE(j,Ll»5,4CFIJ,L»|l*CFIJ,L»2ll 

</24  , 

SSI  I I I -SSI ILI«0T»I9.,ssF( 1 )*l4..SSMU)-5,»SSFIt*l l♦5SFIt♦2l 1/24. 
VSIIII«VSIILI-0T»l9..vsF(|),|9,,ySriul-5,»vSF(L*ll*VSFtL*21l/24, 

4 CONTINUE 
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C*l.L  CO^CT(*,^.^riG,OM)(*,OMXU,l,C,D,cr,N,  r,NP,DrDx,0r0U,W,S,0£U 
print  lOl.  IPI I I Jl , 1*1 ,*) , J, J»1 ,NP,  lO  I 
PRIM  102, (Gill. Mil, W(l|,S(ll,|.l»l.NP.lOl 
1J2  rORHAtlSM  G • I I T7 ,1<(  t I S,  10,2*  I , IS  I I 
in  iLOOP.tOtOlGC  TO  55 
BCTavSl ( 1 l/VSOLO 
55  continue 

VSOLDaVSI I I 1 
DO  5 I»l ,NP1 
5 S(  I I.GI  I l*BET«SI  I I 
DO  2 I I ■ I , 20 

21  DALI  I lal*. I 
)2  DO  8 *■) ,20 

VAL I * I »0. 

DO  * I«l ,NP 

OELllM  aOELI I l-OAL(*l»SI I 1 
* inABSlDELIIIIItGT,OHAxlOEUlII.SIGN(OMAx,OELIIII| 

DO  7 I«1 ,NP| 

I I ■ I * 1 

call  FCT(F,0F0*,X,N,I,dE0U,NP,0EL1I 

DO  RN  Ja I ,N 

RR  X (J, I I la* ( j,  I laOTaF (J,  I I 
CO  TO  7 

RS  do  R5  Ja I ,N 

r5  *IJ,IIIaX(J,II«DTa(55.aE|J,II-59.anj,I-|UJ7,aF(j,I.2I. 

f*.aF( J, 1-31 I/2R. 

call  rCT(F,OFO*,X,N, I 1 ,Of0U,NP,0ELI  I 
DO  Ra  Ja|  ,n 

R4  XIJ,IIIaX(J,II»DTa(f,,F(j,IIU19,,nj,I1.5.anj,I.l|*rtj,I-2n/*Rt 

7 continue 

8 VALUIaEla(x(l,NP).0NaC0S(X(*,NP||-XlFlaa2*TF* 

C E2a|X12,NPU0NaSlNIXI6*NPll-X2F|aa2 

La  1 

VEL I 2 I aVAL I 1 I 
DO  22  Ia2,20 

iriVALI  I I •6E.VCLI2) ICO  Tq  22 

XELl2>aVALI  I I 

Lal 

22  CONTINUE 

IF|VELI2).GT,VAL0IG0  to  23 
IFIL.EG. I ICO  TO  27 
VEL I 1 I aVAL I L- 1 I 
VEL I 3 I aVAL I L* I I 
ALP  I I I aDAL I L- M 
ALPI  2 I aOAL  ID 
ALP  I 3 I aOAL I La  I I 
CO  To  2R 
27  VEL I 1 laVALO 
VEL I 3 I aVAL I 2 I 
ALPI I laO, 

ALPI 2 laOAL I I 1 
ALPI 3 laOAL I 2 I 


«0  TO  »t  i 

|J  ITt»ITE*| 

DO  25  .20 

2S  D*L I I I t I ) /20. 
in  mt'’.»,5)60  TO  25 

60  TO  12  I 

2*  r»INT  I IO|V*LO 

IIO  rOUHATI'  FUNCTION  C0NVER6E0  To  'lEjO.lSI 
60  TO  52 
2T  continue 

RRINT  I02.VAL0, I VEl I II , Ib| ,31 
10)  FORMATITH  VAL  • ,E|5, IO')l%X,E15t lOI I 
IF  I ABS( VALO-vEL I 2 I I ,LT, •001  160  TO  52 

6 I I I I I ■ I • ' 

6 I 2 , I I • I • 

el3>|)>l> 

6( I ,2Ib*LRI  I I 

e(2>2>»*L*’(2l  1 

81 ) t2 I ‘‘LR I 3 I 
81 1 ,3Ib*LPI 1 

8 12,31 "‘LR I 2 I ••2  l 

8<3,3I«*LRI3I»»2 
OIbVEU  1 I-TELI2I 
02«VEU  2 I -VEl  I ) » 

IFIAgSIDl I, LT.. 0000001 1*0  TO  |T 
IFIABSI02I.LT.. 0000001 1*0  TO  IT 
60  TO  20 
IT  ALR0RBALPI2I 
60  TO  I) 

20  call  ROLUIO.OI .VEL,C0CF.3,0,0I 

print  |05,(ALPin,I.I,3l,(COEF(II,I«l,3l 
105  FORMaTITh  alp  • ,3(E15. 10,2X I ,3(E12.7,  1 X I I 
ALPOP»-COEF(?l/(2..C0EFl3»l 
PRINT  lOH.ALPOP, (COEFI I > • I«1 .31 
10"  FORMATITm  alp  • ,5  I E 1 5. > 0 ,5X  1 1 
13  DO  IT  1*1 ,NP 

deli  I IbOELI  I I-ALP0P»SI  I • 
w I I I ■ I • 

IFIABSlDELUII.GT.DMAXtNinBO. 

|T  IF  I AbSIOEL  I I I I .6T.DMAX  iC>ELI  I IbSICNIDhAX.DELI  I I I 

I LOOP* I LOOP* 1 

100  FORMaTIITH  del  ■ »E 15. lO.Al 2X ,E|3.i I ,2X, 15  I I 

101  F0RNAT15H  P • , I T7, 1 1 E iS, 1 0,2x  I , 15  I I 
IFI Il00P.LE.20IG0  TO  16 

52  continue 

print  53, IM,VVAL |P| ,n.] , 101 

53  F0RNAT(10I5H  M • , 15 ,5X »E1 5, 10, / I I 
200  F0RMATII5H  I • ,I6,5X,*HoEL  • ,E15.l0ll 

STOP 
END 
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r 


sue»ouTiNr  rcT t r ,dfdx  ,x .orou ini' lOf  l ) 

OtMtNSION  DrDXlN,NiNF|,*|N,NFl,DFDU(»<,NF|,OEL(><F>,r(NiNPt 

COWMON  Cl,C2,C3|C<*.CS,C*iCr,CS,Cf 

BlaCOSIX(4,MI  I 
B?«'S INI  X I 4,M I ) 

X1»X( I ,M) 

X3»X(2,P) 

X3>X| 3.H) 

XH.X(‘t,M| 

XS-XIS.PI 

l«K«l 

r I I ,H  »«B1 •X3*B2»Xt 
F I 2 iM I •-B2»X  3*B I »XH 
F13,mI»C1»X3«X3«C2»x6»X‘'»C3 
fM,M)«CH»X3»X‘<*CS»X3»x5*C*»X3»X3»OtL<"l 
FI5,Hl»C7»X3»Xt»C6«X3»x5*CB»X3«X3»OEl.lM) 

F ( * 1 M ) »X5 

DFDXI I ,3)H|«BI 

DFDXi 1 ,X,WI.B2 

DrOX|l,4,wi»B2»X3»Bl«X*t 

OFDX(2,3,P)»-B2 

DF0X(  2,",»'I.BI 

DF0X|2|4,"IbB1*X3*B2*X<« 

DFPX(3,3,M|.2t«Cj»X3 
OFOX( 3,t,HI.C2»XS 
OFOX( 3,5,p|«C2»XM 

DFDX|<t|3,P|*C‘)«X<t*CB*XS*2i*C4*X3*DCl-IPI 

BFOX(l,X,F'l»C*t»X3 

0F0X(*t,5,P)»C5»X3 

OFOX(5,3,K)«C7*X*t»CB»XS*2»»C*»XS»OEL(H| 

OFOXtS.I.Pl-CT^XS 
0F0X|6,6,H|.C8»X3 
DF0X(*,6,H|.1 . 

DFDU I B I ■C**X3«X3 
0F0U(5i"I»C9.X3«X3 
RETURN 


ccrcT 


subroutine  COrCT|X,P,PF,GiOHX)(,OHXUiZ,C,OiCriN,M,NR,DrpXfOrDU,W,S, 

*DEL  » 

D1 pension  XIN,NP),P|N,NRliPr(N|NPI|6(NPI,0HXX(N,N|,0HXU(N,NPI 
»,ZI N,NP I ,C  t N.NPI ,D(NP I ,Cf ( N,NP 1 ,orOX(N,N,NP ) ,orOU(N ,NP I ,SINP 1 
l,W|Npl ,DEL(NP) 

common  C!,C2,CJ,CB,CS,C*,C7,CB,C? 

0*T*  INO/0/ 

IMInD.EO.OIPBINT  tO,C|iC2.C3.C<),CS)C*|C7,C8tCf 
10  rORM»T(PH  C • ./f I X ,*1E I 2t 7 , 1 X I I 
P1«P(1|M) 

P2»P|2,M) 

P3.P( 3, Ml 
PN.P I N ,M I 
P5-P(5iM| 

P*»P(t,M| 

XIaXi I •"! 

X2«X I 2 |M I 

X3>X I 3,M I 

X*t«x  ( <t  ,M  1 
X5.X(S,M1 
Xk«X| S,M) 

ei*cosi x*i 

B2»-S1N1X*) 

L-M-1 

PM3,M|..Pl»ei»P2»B2-2,*P3»Cl»X3-PP»(C<t»X‘t*C5»XB»2t»C*»XJ»0EULll 
t .P5*IC7«X<4*C8*XS*2.»C**X3»DELIL)  • 
PF(P,MI«-P|»B2»P2*B1-P3*C2»x5-P‘(«CP»XJ-P5pC7»X3 
PFI5,m)..P3»C2»X<I-PN»CS*X3-PS.CB»x3-P* 
PFIt,MI«-Pl»x3»B2*Pl»X‘l»Bl-P2»X3»Bl»P2»XM»B2 
DHXU(3|M|a;.«x3«|PN«C**PS*C9l 
G(Mla|PN«C«aPS>C«lax3*X^ 

in lND>EQ,OIS(M|aG|M| 
ir(M,tfl.2l IND«2 

0HXX(S»3|a2.«(P3«Cl»(P<(«C**PSaCPI«PELlL)  I 

OMXXlS.PI^Pt^CPaPSaCT 

OMIXiS.^I^PNaCSaPSaCB 

OHXX l3|*)«P|*B2*P2aBl 

OMXX I P , 3 I aOMXX I 3 , P I 

DHXX I P ) S I aP3#C7 

OHXX ( P , 6 I a-P| aBl aP2P82 

OHXX I S ( 3 I aOHXX I 3 t S I 

0HXX|S|P|pDHXXIP,5| 

OHXX |t,3laOHXXI3,AI 
DHXX(B,P|aOHXXIP,4l 

0HXI|t,4la-PialX3*B|*XP*B2UP2P|X3PB2>XPaB|) 

00  I |p3,6 
CE  t 1 ,M1.Q, 

00  2 J>1 |N 

2 Cni,MI»Cril ,M I.OEOXIJ, I ,MlaClJ,MI*0MXX(|,J|P7IJiMt 
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I 


I CM  I ,»)»cri  I ,*»i-OM*u(  I I'SiMi 

0(H).DHXUI3|M)*Z(3,HI*D^0U(<|,Hl*CI'*iM|*DFDUIS,Ml*clS 

RETURN 

END 
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MNlNRL 

COMP  I LEP ( «M«  1 I 

P*R*mETER  NP.lSOl ,NP|»nP-1 ,ML»I50 
DIMENSION  0 I NL I I V I NL  I 

COMMON  PU),«l«i6l,BI*liXltl.GU,*l.lFt6l,*l6,tl|R0l4i6l,X0l6l, 
»U(NP),RSl»,6,NP|,XPL(*,Ni.),i,l|t|,x0lr(6iNP|,XE(6l,X0rSlNP| 

CXLL  COEFCICl.CZ.CJ.CN.Cb.Ct.CT.CS.CPI 
X]p67«S* 

SO  REXD  I , I X I I I , I ■ I ,6  I ,UR,U| ND 
1 FORMAT  I *F 1 O.S I 

REXD  1 , t XF ( I ) , I. 1 ,4 ) , TF 

UP«OP«lOOC* 

p^INT  J3,UIN0 

J3  FORMiTI*  UIND  • ',FiO,5,//| 

print  20,IXI|l,la|,4l,|IFII|,|«|,6l,uP 

20  FORMXTISM  X « |6  ( E 1 5.  10  iSX  I , / ,6m  xF  . , » ( E 1 S . 1 0 , s I I , / , 6M  UP  • , 

•FlOtJ,//) 

DC  1}  l«l ,6 
|3  XlIlKXIII 

X4.XTXN2(XF(J»-X(2I,XF|1)-X(llt 
III,]) mCOS 1X4) 

XI2i3)mS|n(X4) 

XI J.]Im2.«C1*X3 
X ( t ,M I M-X ! 2,31 
X I 2 • <*  I p X I I 1 1 1 
XI<t,l||aC<)PX3 
XIS|P)pC7pX3 
X I •)  ,S  I pCS*X3 
XISiSIpC8*X3 
X I 4 ,S I ■ 1 • 

X I I ,4  I ■-X3«X I 2 . 3 I 
XI2i6IpF3*XI 1,31 
SIRIpC6pI3*X3 
BISIpC«*X3«I3 
I I I I aO, 

I 1 2 1 aO, 

II3laII3l>I3 
II4laXIX|aI4 
TaOt 
lOF laO, 

DO  22  I • I . 2 
SI  I f I iPlOO^ 

RSI  I , I .NPla.iOO. 

22  continue 
1FU«0 

print  2X,CI,C2,C3.C<*,CS(CX,C7,CB,C« 

2*  FORMxTism  C a ,»IE|3. 8,1X1  I 

print  25,  MXM  , JI  ,J-1  ,4'  . 1*1  ,*l  , IBI  I I , |p1  ,41 

2S  FORMXTISM  X a ,*l T4 , 4( t lb. 1 0 ,5X I ,/ 1 ,// ,5m  B a , 4 I C 1 5 . I 0 , 5 X I , / / I 


I P.O 

DO  2 I 1 • 1 , NP  1 
K »NP-  1 1 

♦ 1 

DC  P !■!  |6 
CO  P J» 1 , 6 
“ PI  I iJIpPSI  I , J,p| 

C»l.L  PCEv  I * ,8  ,UP  ,RD  ,P  ) 

DO  1 I ■ I , 6 
DO  } Jp| 

3 PS  I I , J I* I «PS I I , J ," ) -h.bD I I , j I 
IP.IP*! 

ir I Ip.lT.SOIGO  to  2 

3P  CONTIMjE 

PRINT  21,II,»,MPi(I,j,K),j.i;»i,l.i;tl,((B0II,jr,  J.l,41,l»l.*) 

Zl  P0PH»TI4h  II  ■ ,IS,6m  k • , IS , / ,*( * I 5« ,E 1 S. 10  I , / I ,// , 4 I *( S» ,E 1 5f 10 
» I I / I I 
IP»0 

2 continue 

iriUIND.ST.i.lOO  To  22 
00  3l  ««2,NP 
00  3l  !■!  ,4 
00  3 l J«|  ,4 
3 1 PS  I I , J |K I «»S I I , J , 1 I 
32  continue 
IP-C 

00  S 1 •! ,KP1 

1 IpUI 
U I I I aO • 

CO  4 ja  1 , 4 

4 Ul  I laUI  I l•«IJla(PS(a,J,I  laBlPlaPSlSiU,  I laBlSI  I/UP 
00  T ja I , 6 

ID  I J I aB I J I au I M 
DO  7 I a 1 , 4 

7 lO I J I a lO I J I • 4 I J , K I a I ( I I 
UCaU I I I 

DC  1 9 Ja 1 , 4 
|9  PlJIallJIalllJI 

C*LL  ECIP,XE,Cl,C2,C3,C‘',C5,C4.C7.Ce.C9,uCl 
DO  28  ja|  |4 

28  lDIE|J,I|aXElJI>XDljl 

xorsi I |a0. 

DC  29  Ja I ,4 

2*  xors ( I I axOFS  1 I I »X0 I r I J , I 1 aaz 
xor 1 aXOF I aHalOrS I 1 ) 

DO  t Ja|  ,4 

8 XIJIaXIJIaMaxOIJI 

IPalPal 

TaTaM 

IM iPiEQ. lOIGO  TO  9 
60  To  5 
* IPLa|Pta| 
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I 


miNT  IO,U(  I I , I >1  J)  , J«1  (*)  , I 
00  1 i ja  I < * 

II  «PLIJiI'U»*IJI»XllJl*(»ftJI-»I(J)l«T/Tr 
POINT  1 8 , T , I » PL ( J , IPl  I , J«1  ,6  I 
)8  rOON^TIBM  time  p .t 1 O.S •*( 2x ,E is. 1 0 I I 
IPpO 

lO  rORM*Tt*M  U P ,E I S. 10, J* .* I E 1 5. 10. I « I . IS  I 

5 continue 

U I NTaO , 

00  27  iPl ,NPl 
27  uINT.UINTpOMIpn 

PRINT  30,UlNT,J0rI,|I,»DfSIM,IX0IPIj.II,Jp|,6l,Ipr.NPl,10l 
jO  rORM«TI8M  UINT  p ,E15.|C,8M  XDPI  p ,E15.lO,/,8M  XCIP  p ,/,I|X,I8, 
»IX,7(E1S. 10,2X1  I I 
IPpO 
IPLPO 

00  12  iPl.NL 
01 1 IpXPLt  1 , I I 
I 2 VII  ) pXPL 12.11 
C*LL  CPNPlT 
CALL  XTEPlT 
CALL  LINPlT I 0 , V ,NL ) 

00  lA  iPl  ,A 
|R  XIIIpXIIII 

00  15  iPl ,NP| 

UCpU I I I 

CALL  PCIX,X0,Cl,C2,C3.C‘',C5.C*,C7,Ce.Cf.uCl 

00  U JPl  ,A 

I A X IJ I pX IJ I pHpxO I J I 
IPpIPpI 

in IP.EO. 10160  TO  17 
50  TO  IS 
l7  1PLPIPL»1 

print  1 0,U I n , I X I j I , j, 1 , t I , I 

01  IPlIp*'  1 ' 

VI IPlIpXI2I 
IPpC 

l5  continue 

call  CU»VEIO,V,Nl,OI 
call  enopl  1 1 ) 

STOP 

ENO 
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COEFC 


SfBRCUTlNE  C0EFC(Ci,C2,C3,Ch,C5,C4,C7.C8,C9I 

RH.2S* 

RH0»2 
RL« I S . 

XT»..7,S 
RiC-l  .07/2.  1 
DP».  1 5 

RF.R|,/  I 1 • 1 7e»SS>»T  I At  I I 

RK1*.65«RF»»i-|.m9) 

BA2*.73«Rr»»(.ll3'>l 
R>*L"8  00.»I  l.»RA|  )/J2.2 
R'‘T*eOO.»(  1 .♦Rk2)/32.2 
RrR«.<*2»RF«».331 
RJ2«ROO.*3RO.«Rap 
RLTAP.2.S 

THRUST»(RH0/2.)»AA.r*,47,54,(,7,^4 

m8s(RH0/2.  )«AA*R|,TtP 
m2> I RHO/2. I *AA 
H 1 •h2»DR 
h3»"2»*L» • 3 
1 . 3»M2»RL 
m7«-,7»H2»RL«RL 
C1»-m1/Rmu 
C2»R>*T/Rml 
C3«ThRuST/RML 
C*-  ( m2»m8  t /R*«T 
C8pIm3*P*U.'*’8»»Ta|/pmt 
CA«Hg«RK /PMT 
C7» I m*-XTA»H8 I /R J2 
C8«(M7-XTA»XTA«Hgl/Rj7 
C9»XTA«Mg.RA/Rj2 
return 

ENO 
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I 


rc 


suBBouTii.r  f c ( « , »D , c 1 ,c? . 

DIxCmSION  X(6) •X0l6) 

XDinaX<]l>C0S(>l6ll>l('*l 

XCl2l»X(3l»SlMX(*)|,X(») 

XO|3|»C1»XI3I«»?*C2»XImI» 

x0|‘*I"C«»X(3I»X(<41»CS»X<3 

xD(5l»C7»X(3l»xii(l»Ce*x'3 

I 0 ( « I • X I S I 
« C T U P N 
END 


C3,C‘t 

.CE,, 

C6  , C7 

.CS 

,C9 ,u ) 

•S  1 N ( 

X ( 6 ) 

1 

•COS  ( 

k 1 6 ) 

1 

X 1 S 1 • 

C3 

1 »x  (s 

1 «C6 

•k(3I 

•X  ( 

3 1 »U 

1 •!  1 S 

) ^C9 

•X  ( 3 ) 

•X  ( 

3 1 *U 
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fCtv 


SUB' OUT  I NE  E ct  V ( 

> . B , UP , "c , R 1 

0 I wEiiiS  1 ON  *(  « ,4  1 

,B  ( 6 1 ,r0 ( 6 ,6 1 

00  1 

J • 1 1 6 

00  1 

J*  1 1 * 

eo  ( 1 

•6(R|.rI|,si« 

00  1 

>•1.6 

*0  ( 1 

, J 1 •«P 1 ! , J 1 - 

* 1 » . I 1 •« l> , j 1 

Ot^L'SN 

ESC 


knlnep 


I 


i 


CC«P I lE  ® I >"•  1 I 

double  precision  HGRwo.NaPwi 
OOUBlE  precision  D,E,E2 
DOuBlE  precision  COD.CO^.CChI ,cood 
OOLBlE  precision  mgrw.hDmC.hDpO 
• *R»>«ETER  N0P«I3 

p1R*hETER  N«6,*‘«|2,NP.5I,I0I"N»NP*IiI02«N*>»*S,N2*n-2 
OInEnSIOn  El  1021  ,E2I  102* 

OIPEnSION  XPLTINP) .TPlTInP) .XPTgiNP) ^TPEgINPI 
C0M“0N/X/0rDCX(4,6,6l  iD'JTEST 

C 0 X « 0 N / C X L L / T r 

COXNON  XINI  .XlNINI  ,XPL|R,I80I  ,*0INI  ,*0(Nl  .XOOINI  ,OXD*(N'm)'f(N)  , 
loxoox  IN,X)  .OxDXOlN)  ,0»Dt>TlNP,Nl  ,X|N,m)  .CEUNP  I .OPlDll  I lOI  , 1021  , 
fOxDDiCINI  ,0XDT|N) iDxDOTIM  ,oro>IN,NI  ,N|Np,NI  ,W1  I lOI  I ,*1  I 1021  , 
»D»D*|NP,N,N,h| ,DnOxO(NPiN,8)  ,m6RnI  102,1021  .mDxci  IDD  ,mDmO(  ID2I  , 

»mGR"III02, 102), OtIC2l,D'’OT|NI,XFP|N),Xor(N|,X|0(ID2), *151102, 301. 

»*1PII02.2B),CGC*(I02),C'’CIN, 180), C0B(X, 1801, COO(M),COH|m,M|, 
fCOxl |X,"1 .CODCIP) , COGIN, M) 

IN*.  1 
XFP I 1 ) .0 . 

IFP I J I .0 . 

PI.3, lN15»2i536 
COEF,  1 , 

ON«e 

22  FORMtT I 8F 1 0 , M I 

print  23 ,C0EF ,XT0 ,YT0 ,*Lph* 

23  F0RM*TI8H  COEF  . ,F10,R,6HXT0  • , F I 0 . H , 2 X , 6 H T 0 « , F 1 0 . R , 2 X , F 1 0 . 2 I 
ZETEsT.SGRT I 3.1/2. 

COEEO.COEF 
0UTEST..0C02 
Dn*X.15./S2,2»S76 
CUTES2.Dm»x 
21  CALL  C*LLIN 
DO  13  I«l  ,N 
F I I 1.0. 

0X0*0  1 1 1.0. 

0X00*0  1 I 1.0. 

0 X OT ( I 1 «0 . 

OXOCT I I I *0. 

DO  I R J. I , N 
0X0*1  I ,J).0. 

0X00*1  I , J ) >0 . 

IR  C06I I ,JI«0, 

X 0 I I I .0 . 

*0F I I I .0. 

XFPI  I 1.0. 

00  2S  «"1 ,N 
2B  OFOX ( I ,r 1.0. 

1 3 orOT ( I 1 .0. 
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I 


00  IS  1*1." 
rOD  M I "0, 

coco ( I ) "0 . 
on  IS  j* 1 1 « 

|S  COHl ( I , Jl*0. 

00  U t*l 
CEU  I I *0  . 

CO  16  J • 1 1 ^ 

DWCDT ( I , J I *0 . 

* I I I J I *0 ( 

00  17  K»1  ,8 

|7  C»D»0II,JiK!.0. 

00  16  » • I . *• 

00  16  L«  1 , 

1»  OfcO*(l,J,i,LI*C. 

00  18  1*1 . 102 
CG06 I I I *0 • 

All  I 02  I *0 . 
hOfC I I 02  I «0. 

HOMO  I I 02  I . 

01  I021-0, 

00  2h  J* 1 , 28 
A I S I I I J I >0 . 

A1P(I,JI»C, 

A1SII,2’’I»0. 

A 1 S I t I 30  t.O. 

A 1 0 I I I «0 . 

00  18  J* I t I 02 
hG8p(  I |J1*0» 

I • mGRx I ( I I J I «0 , 

00  19  I • 1 , 10  1 

» 1 I 1 0 I I «0 • 

00  19  J* 1 , 1 02 

1»  OxlCAIll.JI.Ot 
DT.Tf/ 1 79. 

00  IB  I • I , 6 

AOIl|«XIMI) 

9«  Aconi**iii-»i*ii[) 

00  12  «■!  , I 80 
T« ( A - I . I *07 
00  JA 1 ,1 

93  C0CIJ,Al9»PLIJ.AI-A0IJ|-AD0IJI»T/Tr 
00  9«  J9  1 I H 

99  COB  I J 1 A ' 9^ I N 1 J9^ 1 »T /7F  ) 

92  CONTINUE 

00  HC,  I9J  ,M 

00  95  J»1 ,M 

COh I I , U I «Q  . 

00  95  *91,180 

9S  COmII,JI.COhiI,ji*COBIIiA)9COBIj,kI 
CALL  pOLUICOm,COhI,COO,OoOC,n,C,  1) 
DO  96  I 9 1 . N 


DO  •<*  J«  1 , H 
COGI I ,JI*0, 

DO  It  . I , 1 8C 

“6  COCII,JI’COGIt,JfCOCIIiKltCOpij.Kl 
DO  «7  I ■ I >N 
00  <17  J«  1 , M 
t ( I , J I aOt 
00  tT  K* 1 ,H 

h7  k I i,ji»coGi  I ,iti»coHnK,ji 

eRINT  |05,TF,(*0(I),I.lt6l,(*00(Il.l.l.5) 

OT.Tr/  I ►tP*  I t I 
DKG..  1 
FRS-.Ol 
E RS« . 00 1 
ILOR.O 

■lO  I lor.  I1.0P*  1 
I lOOp.O 
OUG.DItG/ID. 

PRINT  MliEPSiDUG 

Nl  r0RPjT|7H  EPS  ■ ,020. ISib* ,4NDI » . ,E20,l5l 
170  fcr-«t ( trro. 1 0 1 

PRINT  l03,t(»ll,JI,J.r.P),I.l,N) 

20  SGDE.SOPTtOT/EPSl 
DO  SOO  I . 1 , I D2 
00  ScO  J» 1 , I 02 
SOD  hGRn I I , J I pO . 

*T*Rg»1000. 

T T »Rg»0  « 

DO  I L»1 ,NP 

C*H-  0CT(»,io,*,»0,»OO,OxC*,OrDD*,DX0A0,D«DO»O,DXDT,OXDDT,L,N,M,NP 
l,DT,iLPNX,XOF  I 

call  OPCT(r,*,XD,OrDX,L,N,OT,DEL,NP,iTD,TTOI 
XPLT(LI*«I 1 I.XTARG 
TPLTiL l■X(2)♦TTARG 

xPTGILIpxtabg 
TPTGiL'ptTARG 
90  CONT I nuE 
DO  2 J*  1 ,*■ 

DO  2 X . I , N 
DO  2 I . I , N 

0»DA|L,J,',l|p-Dr0X(j,Kl«CX0A(X,I)»S0CE 
2 IF(J,E0iKIDpDA(L,J,x,I|P0N0AIL,J,X,I|9CXDDA(J,II.SQ0E 
CO  65  jp 1 , N 
00  66  >P|,N2 

DP0A0iL,J,Xlp-Dr0X(j,r«2|.0x04ClK«2I.SaDE 
66  ir(J,E8.X‘2IOPDAo(L,J,X>pOpOAOIL,J,XI»SSDEpOXDOAO(JI 
65  continue 

DO  87  JP I ,N 
DEDT ( J I pO. 

DO  87  1 P 1 ,N 

8T  DFDT(J|pDFOT(JI»DECx(J,I  IpOxOTI  I I 
00  3 Jp  1 ,N 


1 

1 
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1 


1 


m 


I 


kil  . j)  •'  xtH  j I ijt ) .sacf 

3 0*0t’TIL,JI**ILiJ)/(2.»D'f|»(0iC'DTIj)-0f0T(Ji  l•S(;OE 
TIMt.(L-l . (•ST 
50  To  ^3’ 

in  IlOO^.LT .5  I GO  To  50» 

00  501  I?*l tN 
00  501  I***!  >** 

j I • I I H 

00  50?  I3«l 
00  50?  15*1  •** 

J?«  I I 6 

00  502  I 1 • 1 i ^ 

502  MGSfciJl ,J2I*m5®«(Ji  ,j2)-0E0DXI  1 j ,12,13I»D«D»I  I2,Im>»CiD*I  I3,I5I« 
f S:OE  I L , I 1 I 
00  503  I3»l ,« 

I5«N«^* 1 3 
00  505  I 1 • 1 

505  m&BwIJI  ,J3)«mG®W(J1  ,J3I*DED0xi  I [ , 12, I31»0X0*I  I2,I«I»0X5»0(  I3l«S00E 


f •« I L , I M 

50  3 mGOi*  I J3  , J I I •HG'"  ( J I , J3  1 
DC  So**  I 1 • 1 ,N 

50**  M6U>ilJI,IO2l*HGBt<IJl,|02|»(0*i0*(L,tl,l2,l‘<)/l2,»0TI-OX0O»ll2,I‘'l* 

iS5DE/0T  I .k ( L , I 1 I 
501  mG»» I I 02 , J 1 I .hGK w I J 1 , I 02 ) 

DO  50*  J2« 1 , H 
j3"^«’^*02 
00  5CT  I 1»1 ,N 

50?  xG*»(J3,IC2I»hG»w(jj,Id2)»(0*<CX0(L,Ii, 121/12, •OT|.OXOO*0<J2»21» 

f SGOE/DT  I •!»  I L , I I I 
50*  «6">' I I 0 2 , J3  I .HGSw  ( J3  , 1 o2  ) 

DO  5CB  I I • 1 , N 

50*  mGBui  ID2,ID2)»*<G(!i.(  ID2,IC2I*(S3DE»X0(  I i l-wM.,!!  l/XtMOIL,!!  I/DT/DT 
50*  CCKTinuE 

POINT  1 05, time, (x(,), 1. 1, 4,^,, cm, 1.1,61, mil, I, i;*i 
105  F0»h*T(5h  T • ,riO,6,HHX  ■ , 6E 1 5 . 5 , / , *x  XO  • ,6E2c,5,/,5m  r • , 

< *E20<S> 

» T X 0 G ■ * T * 0 5 ♦ « T 0 • 0 T 
»T»Bg*TTXPG»TTC«DT 

I continue 

P»INT  50,11  ,0EL ( I 1 , I • 1 ,Np I 
50  rOPP*Tl5N  I . , I3,2X,6mD£i.  . ,E20.15) 

G.SQpT ( ( nP-  1 . I .DT I 
OGOT.G/ t 2. xOT I 
DO  •*  J.  1 , N 
DO  * I«1 ,« 

<<  111  J.N-x.  I I .6  I J , I I 
*111021 *0T 
00  *7  !•! ,N2 
*1 IN.H.t l.*00l 1*2) 

67  CONTINUE 
COEE.COErO 
00  IN*  1.1,102 
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r 


I <|«  < 1 0 I M 1 I I ) 

00  5 J. I , N 
DO  5 L«1 

5x11 J.NP-nP*L 1 ■» I L I J ) 

XI  I IDI  !■& 

DO  « Jxl  ,S 
DO  « Lx  I ,NP 
DO  7 Kxt  ,N 
DO  7 Ixl  ,M 

7 DX1D*1IJ«NP.nP»L,k«M-M»I|x0X0a(l,J.K.II 
DO  68  Kx| ,N2 

*8  OWIDaI (J«NP-NP»L,N»M*A)*0PDA0IL,J,AI 
4 ow  1 Da  I ( J«NP-NP»L  , I D2  I •DWODT  ( i, , j | 

DK I DA  1 I 1 D I , 1 D; ) xDSDT 
DO  9 Ixl , ID2 
DO  8 Jx 1 , I D2 
DC  8 Kxl , IDl 

8 M6PX(  I (JlxHCPXI  I ,JV»0xi0a1  IK,  I IxDXIDai IK ,JI 
MO«C I 1 I xQ, 

DO  9 Kx 1 , I D 1 

* mOHCi  I IxxDMCl  I I-OWIOAI |K,  I IxW]  IK) 

CVxO, 

DO  ir  Ixl , IDI 
I 2 CVxCvxwl  I II X.2 

1 TxO 

c voxcv 

print  102, iloop.cvo 

)D2  F0P«aT|5M  I X , 16 , 2k ,EjO, 1 5 I 
DO  3x  I X I , ID2 

JX  MGRx I I , I 1 xMGPx I I , I 1 »0 I a6 
GO  To  xOJ 
DO  X02  1x1,102 
DO  XQZ  Jxl, 102 
«02  «GPW I I 1 , J I xMGPx I I , J 1 

call  TRED 1 I ID2 . I 02 ,mgrx I ,0  ,E  ,t2  1 

C*LL  ImTgl  I M 02 ,0,E , IErR I 
PRIM  xOO  , lERR  , I D I I 1 , I . I , ID2  I 
XOC  rORPATI7HlEP  X ,I5,xho  X , 1 00  I T 1 8 ,E2o.  1 5 ,/ I I 
X03  continue 

C*LL  RDLUIHGRN,hgRxI ,H0"C ,HOHO, ID2  .0.01 

GO  To  7* 

76  DO  7h  1x1,102 
7X  MONO  I 1 I XHPHC I I ) 

74  continue 

00  10  1x1,102 
4 1 I n xA 1 I I I xCOf PxMOMO I I ' 

I 0 CONTINUE 

DO  11  Jxi  ,N 
DO  II  |xi  ,N 

II  AU,  tlxA|  luXH'Hxt  I 
00  4R  |xl ,N2 

ADO  I I >2  I xA I I NxH* I I 
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I 


4*  continue 

DT«* 1 (102) 
ir (OT.U^.O.  I0T«.01 
print  1 03 , I ( 4 ( I . J I , J. 1 , H I , I . I ,N I 
lOJ  rORN*TI*)IH  ,*E20. I5,/,T3i6E20. I5|/I .//) 

cv»o; 

00  I S I L- 1 ,NP 
S00F.S5»T lOT/EPS) 

C»LL  ECT(*,XD,»,*0,A0O,0*D*,D»0O*,O«0t0.D«OO*O,0X0T,O)(00T,L,N, 
I , OT , alpha , AOF ) 

call  OrCT(r,x,*0,OFD«,LtN,DT,DEL,NP,xTO.TTD) 

00  ISI  U*1 >N 

* ( L I j I • ( »o ( j ) -r ( j ) I •sapt 
CV»Cv«w(L,JI»»? 

151  continue 

6»S0PT I ) NP. 1 , ) »DT ) 

CV*Cv*G»G 
GO  To  153 

IE(Cv*LF.CV01GC  to  i53 
COEF.COEF/m, 

00  158  I • 1 I 1 02 
158  A 1 ( I I >A I 0 ( 1 I 
I T« I T* 1 

IFI 1T.GT,5)G0  to  153 
GO  To  7* 

153  continue 

print  ioh.iaoi 1 > ,i«i  ,n) 

PRINT  1 OR , ( A00(  I ) , I. 1 ,N I 
1 OR  FORM  at ) 1 M , 4E20 . 1 5 I 
DE"Cv-&*G 
TF* ( nP- I . I »0T 
I LOOpr I loop* 1 
PRINT  100, ILOOP.TF.CV.D^ 

100  F0Rh*TI5h  1 ■ ,Ir,2x,5m7f  • , E I 5 . 5 , 2 X , 5 hc V ■ , E 1 5 . 5 , 2X , 5HDE  ■ 
»E 15.5) 

00  UO  1>1  , 102 
140  A 1 S ( I , I LOOP ) rA 1 ( I I 

IF ( AgS ( C VO-CV I .LT . ,000 1 • GO  TO  R20 
IFI IlOOP.LT. ISIGO  TO  20 
R2D  IL0P2RIL00P-2 

145  F0Rh»T17h  41P  R ,(T8.Ei5,10M 
TPPr.OI 
TrO. 

DTr.oOI 
00  120  Ir1,N 
120  XII  Iraqi  i i 
LrNP* 1 
LIhrTF/.OOI 
XPTGI  1 |RX  I 1 l*XTARG 
TPTGl 1 |PX(2I*TTAR6 
TTESTrTF/INP-1  . 1 
TTESTOrTTEST 
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I 


I 23,T, ( XI J) , j.j ,N> 

00  121  I»1  ,LI*< 

If  ( T.  CT.  TTEST  IL»l.»  1 

If (T ,6T,TTESTIXPTGIL-MP)»X| j )*XT4PG 

IflT,GTtTTEST)TPTG(|.-NPI»X(2l'»>'T»R6 

If|T,fiT.TTEST|TTEST,TTEST*TTESfO 

CXUL  OfCT(f,x,XO,OfDX,l.iNiDT,OEL,Nf|XTO.rTO> 

00  122  JPt>N 

122  IIJIaXIJI*OT*f(JI 
T.T»0T 

If IT.LT.TPPIGO  TO  121 
TPP»TPP»,01 

print  1 23,T  , I X IJ)  , j. j ) 

123  fORPxTIGH  T . ,f 10.3, *ElS. 101 
121  continue 

XPT6 ( NP I »X I I ) *XTxRs 
YPTG(NP)aX(2l«YTiRG 
PRINT  I 23 ,T , I X I J) , j. I , 4 ) 
c*Lt  opnplt 

C*LL  XYEPLT 

CALL  L INPlT ( XPLT , YPlT ,nR| 
call  curve  I XPTG , YPTG.NPtO I 
call  ENDPL I 1 1 
91  continue 
GO  TO  21 
230  continue 
STOP 
ENO 


C*LLIN 


COMP  I LtP  ( XM« I ) 
subroutine  C*LLIN 

P*R*»ETER  NP.IBOI.NPI.nP.J .NLpIBO 

COMMON/C*LL/Tr 

0 1 mCnSI ON  Ot  NU ' , V I NL ) 

COMMON  «r(6l,»M4l,XPLI‘,NL),H6,*)tBU',Xl6),G<*i*l,RI*»,Rl*,*l, 
$RDI4,6),«rM*),UINP|,PS(*,t,NPt,X0Ift*|NP),XEU),X0FSINP| 

C*UL  C0EFCICl,C2,C3,C<t,Cs,C»,C7,ce.C«l 

X3P67.S6  - 

50  READ  1 , ( X I n , I p I , 4 I ,UP ,Ul ND 
1 F0RM*T18F10.5) 

READ  1 ,tXF(  I)  ,1*1  ,4)  ,TF 
UPpUPp 1000. 

PRINT  3J,UIND 

33  F0RM4TI'  UINO  p ',F|0.5,//I 

print  20,(X(II.Ip1,4I,(XF1II,Ip|,4I.uF 
20  F0RMXTI5H  X p ,4  I E 1 5. 1 0 *SX  I , / , 6H  xF  ■ , 4 ( E I 5 . 1 0 . 5 X I . / , 4M  UP  p , 

»F 10. 3 .//  1 
DO  13  I«1  .4 
13  xl ( I |pX( M 

X4pATAn2IXF(2I*X(2I,XF(1|.X(1)I 
A( 1 ,3IpC0SIX4> 

412,3) PSINIX4I 
A(3i3>p3.p01pX3 
All, N )P-A I 2,3) 

AI2,h)pA| 1 ,3) 

AIM,i«IpC<*pX3 
ais,m)pC7pX3 
4 I P ,5  I pCSpX3 
AI5,5)pCBpX3 
A I 4 , S I P I . 

Al I ,4>p*X3pA|2,3I 
Al 2,4)pX3pAI 1,3) 

BIp)pC4pi3pX3 

BI5IpC9>X3pX3 

III ) pO. 

X I 2 I pO. 

I 1 3 ) pX I 3 ) >X3 

I I 4 ) pX I 4 ) -X4 
TpQ. 

XOF I.O. 

00  22  !■!  ,2 

ei I , I )pioo. 

RSI  I , I ,NP|P.100. 

22  continue 
IRI-PO 

print  24,Cl,Cr,C3,CR,C5iC4,CT,C8,CR 
24  F0PM4T(5m  C p ,»IEi3,»,1xI  I 

PRINT  25,MA(|,J|,j,|,*),|,,,*|,,gll,,i,,,4| 


25  FORMiTISM  i ■ I * i t * , 1 1 [ I 5i  1 0 ,S)t  I i / I I// , 5m  B • , * ( E I 5 . 1 C , 5 » l',  / / I 
m.TF/(NP.|  . I 
|PaO 

DO  2 I I»l ,NP| 

K«NP.1J 

M»K  ♦ 1 

DO  P I-l  .* 

DO  <•  J*  I ,« 

" PI  I 1 jl»"SI I , J,P) 

C*LU  FCE*(»,B,UP,»D,P| 

DO  3 l-l .* 

DO  3 I I 6 

3 RSI  I , J ,K  l-PSt  I I J,"  )-M.pjO|  I , J I 
|P«IP«1 

IFI IP.LT.SOICO  To  2 
3*1  COsTtNUE 

•RIM  21,II,lt,IIRS|l,j,»),j,r,  t|,!.i;ft|,((RDII,Ji;j.l,6l,I-l,»l 

21  F0RM»T(6M  II  ■ ,I5,5H  « ■ , 16 , / , 61 ‘I 5X ,E iS. 10  ) ,/ ( .// ,6  I 6 I 5X ,E 1 5.  1 0 

S  I • / I I 

IP»0 

2 continue 

IF  I UIND.GT. 1 , I SO  To  32 
DO  3l  F«2,NP 
DC  3|  1*1  .* 

DC  3 I ja I I 6 

31  RSI  I ,J|FI«RSl I ,J,t  I 

32  continue 
IP«0 

DO  5 I»l ,NP1 
I ImUI 
Ul I )>0> 

DO  6 U« I I 6 

6 UIII»UIII»*IJI»IRSl<(,j,I|«B(<t)*RSi5,j,II»B(5l)/UP 
DO  7 JM I , 6 

XD I J I aB I J I *U I I I 
DO  7 Fa  1 , 6 

7 XDI UlaXOl J l»*IJ,X)axlF) 

UCaU(  I 1 

00  1 R Ja  1 , 6 
I*  PlJIallUlaXIIJl 

C»LL  TCIP.XE, Cl, C2, 03,0*1,05.06,07. 08, CR.UC) 

DC  28  Ja|  ,6 

21  XOIF(J,I)aXE(JI»XD(JI 
XDFS I II aO. 

DO  2R  ja),6 

2*  XDFSl  I laxOFSl  I laxOlFI J, I laa2 
XDFIalDF|aMaxOFS(  I I 
DO  8 ja I , 6 

8 X I J I aX I U I aHaiOl J I 

IPalPal 

TaTan 

IF  I IP. ES. I 0 I GO  TO  f 
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60  TO  5 

print  I 0 ,UI 1 ) I I X IJ I , ja I ,4  I , 1 
DO  II  ja I ,4 

II  »PLlJ,|Rl.I"«(JI**I(JI*l«F(J)-«I(j||aT/Tr 
print  I 8 ,T , ( xPL ( J , IPL 1 , Ja 1 ,4  I 
|»  rORH»T(8H  TIPE  a i £ 1 0 . 5 t 4 I 2 X , t I 5 . I 0 I I 

IPaO 

|0  rOR"*TI4H  U a 1 1 I 5 . I 0 , 3 * i * I E I 5 . I 0 i I « I i I S I 

s continue 

UiNTaOa 
00  27  I • I ,NP 1 
2T  U I NTaU I NT»U ( I I »M 

PRINT  30,UINT,X0rl,II,x0FS(II,IX0IFIJ,|l,Ja|,4l,laliNP|,|01 

SO  fORN»T(«N  UInT  a ,E|5,iO,#h  XOEI  a iEI5.|0./,RH  XOIE  a ,/,(|X,It, 
t|X,7|E|S. 10, 2X111 
IPaO 
IPL-0 

00  12  I ■ I ,NL 
0 I I I aXPL I I • I I 
I 7 V I I I aXPL I 2 , I I 
C*LL  opnplt 
CALL  XTEPLT 

call  l inplt 1 o,v  ,nl  I 

00  IH  !■! ,A 
I « X I I I aX I I I I 

DO  IS  l•l,NP| 

UCaU I I I 

CALL  fClX.XO.CI  ,C2,CS,CR,CS.C4.C7.C8.CR,UC| 

00  I4  ja|  ,4 

I A X t U I aX I J I *H>X0I J I 
iPalPal 

in  I p>ca.  1 0 1 GO  TO  17 
GO  TO  16 
|7  IPL«IPL*I 

PRINT  1 0 ,Ut  I I , I X ( J I , Ja I ,4  I , I 

01  IPl I •* I I I 
VI IPlI«XI2I 

IPaO 

l5  continue 

call  CuRvE ( 0 ,V ,NL ,0 I 
call  ENOPL  I I I 

RETURN 

END 
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